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: G4BraggIonModel.cc,v 1.22 2008/10/22 16:00:57 vnivanch Exp $ >> 27 // GEANT4 tag $Name: geant4-09-02-patch-03 $ >> 28 // 26 // ------------------------------------------- 29 // ------------------------------------------------------------------- 27 // 30 // 28 // GEANT4 Class file 31 // GEANT4 Class file 29 // 32 // 30 // 33 // 31 // File name: G4BraggIonModel 34 // File name: G4BraggIonModel 32 // 35 // 33 // Author: Vladimir Ivanchenko 36 // Author: Vladimir Ivanchenko 34 // 37 // 35 // Creation date: 13.10.2004 38 // Creation date: 13.10.2004 36 // 39 // 37 // Modifications: 40 // Modifications: 38 // 11-05-05 Major optimisation of internal int 41 // 11-05-05 Major optimisation of internal interfaces (V.Ivantchenko) 39 // 29-11-05 Do not use G4Alpha class (V.Ivantc 42 // 29-11-05 Do not use G4Alpha class (V.Ivantchenko) 40 // 15-02-06 ComputeCrossSectionPerElectron, Co 43 // 15-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma) 41 // 25-04-06 Add stopping data from ASTAR (V.Iv 44 // 25-04-06 Add stopping data from ASTAR (V.Ivanchenko) 42 // 23-10-06 Reduce lowestKinEnergy to 0.25 keV 45 // 23-10-06 Reduce lowestKinEnergy to 0.25 keV (V.Ivanchenko) 43 // 12-08-08 Added methods GetParticleCharge, G 46 // 12-08-08 Added methods GetParticleCharge, GetChargeSquareRatio, 44 // CorrectionsAlongStep needed for io 47 // CorrectionsAlongStep needed for ions(V.Ivanchenko) 45 // 48 // 46 49 47 // Class Description: 50 // Class Description: 48 // 51 // 49 // Implementation of energy loss and delta-ele 52 // Implementation of energy loss and delta-electron production by 50 // slow charged heavy particles 53 // slow charged heavy particles 51 54 52 // ------------------------------------------- 55 // ------------------------------------------------------------------- 53 // 56 // 54 57 >> 58 55 //....oooOO0OOooo........oooOO0OOooo........oo 59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 56 //....oooOO0OOooo........oooOO0OOooo........oo 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 57 61 58 #include "G4BraggIonModel.hh" 62 #include "G4BraggIonModel.hh" 59 #include "G4PhysicalConstants.hh" << 60 #include "G4SystemOfUnits.hh" << 61 #include "Randomize.hh" 63 #include "Randomize.hh" 62 #include "G4Electron.hh" 64 #include "G4Electron.hh" 63 #include "G4ParticleChangeForLoss.hh" 65 #include "G4ParticleChangeForLoss.hh" >> 66 #include "G4LossTableManager.hh" 64 #include "G4EmCorrections.hh" 67 #include "G4EmCorrections.hh" 65 #include "G4DeltaAngle.hh" << 66 #include "G4ICRU90StoppingData.hh" << 67 #include "G4ASTARStopping.hh" << 68 #include "G4PSTARStopping.hh" << 69 #include "G4NistManager.hh" << 70 #include "G4Log.hh" << 71 #include "G4Exp.hh" << 72 #include "G4AutoLock.hh" << 73 68 74 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 75 70 76 G4ASTARStopping* G4BraggIonModel::fASTAR = nul << 71 using namespace std; 77 << 78 namespace << 79 { << 80 G4Mutex alphaMutex = G4MUTEX_INITIALIZER; << 81 } << 82 72 83 G4BraggIonModel::G4BraggIonModel(const G4Parti 73 G4BraggIonModel::G4BraggIonModel(const G4ParticleDefinition* p, 84 const G4Strin 74 const G4String& nam) 85 : G4BraggModel(p, nam) << 75 : G4VEmModel(nam), 86 { << 76 corr(0), 87 HeMass = 3.727417*CLHEP::GeV; << 77 particle(0), 88 massFactor = 1000.*CLHEP::amu_c2/HeMass; << 78 fParticleChange(0), >> 79 iMolecula(0), >> 80 isIon(false), >> 81 isInitialised(false) >> 82 { >> 83 if(p) SetParticle(p); >> 84 SetHighEnergyLimit(2.0*MeV); >> 85 >> 86 HeMass = 3.727417*GeV; >> 87 rateMassHe2p = HeMass/proton_mass_c2; >> 88 lowestKinEnergy = 1.0*keV/rateMassHe2p; >> 89 massFactor = 1000.*amu_c2/HeMass; >> 90 theZieglerFactor = eV*cm2*1.0e-15; >> 91 theElectron = G4Electron::Electron(); 89 } 92 } 90 93 91 //....oooOO0OOooo........oooOO0OOooo........oo 94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 92 95 93 G4BraggIonModel::~G4BraggIonModel() 96 G4BraggIonModel::~G4BraggIonModel() >> 97 {} >> 98 >> 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 100 >> 101 G4double G4BraggIonModel::MinEnergyCut(const G4ParticleDefinition*, >> 102 const G4MaterialCutsCouple* couple) 94 { 103 { 95 if(isFirstAlpha) { << 104 return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy(); 96 delete fASTAR; << 97 fASTAR = nullptr; << 98 } << 99 } 105 } 100 106 101 //....oooOO0OOooo........oooOO0OOooo........oo 107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 102 108 103 void G4BraggIonModel::Initialise(const G4Parti 109 void G4BraggIonModel::Initialise(const G4ParticleDefinition* p, 104 const G4DataV << 110 const G4DataVector&) 105 { 111 { 106 G4BraggModel::Initialise(p, ref); << 112 if(p != particle) SetParticle(p); 107 const G4String& pname = particle->GetParticl << 113 108 if(pname == "alpha") { isAlpha = true; } << 114 corrFactor = chargeSquare; 109 if(isAlpha && fASTAR == nullptr) { << 115 110 G4AutoLock l(&alphaMutex); << 116 if(!isInitialised) { 111 if(fASTAR == nullptr) { << 117 isInitialised = true; 112 isFirstAlpha = true; << 118 113 fASTAR = new G4ASTARStopping(); << 119 G4String pname = particle->GetParticleName(); >> 120 if(particle->GetParticleType() == "nucleus" && >> 121 pname != "deuteron" && pname != "triton") isIon = true; >> 122 >> 123 corr = G4LossTableManager::Instance()->EmCorrections(); >> 124 >> 125 if(!fParticleChange) { >> 126 if(pParticleChange) { >> 127 fParticleChange = >> 128 reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange); >> 129 } else { >> 130 fParticleChange = new G4ParticleChangeForLoss(); >> 131 } 114 } 132 } 115 l.unlock(); << 116 } << 117 if(isFirstAlpha) { << 118 fASTAR->Initialise(); << 119 } 133 } 120 } 134 } 121 135 122 << 123 //....oooOO0OOooo........oooOO0OOooo........oo 136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 124 137 125 G4double G4BraggIonModel::GetChargeSquareRatio 138 G4double G4BraggIonModel::GetChargeSquareRatio(const G4ParticleDefinition* p, 126 << 139 const G4Material* mat, 127 << 140 G4double kineticEnergy) >> 141 { >> 142 // this method is called only for ions >> 143 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy); >> 144 corrFactor = q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy); >> 145 return corrFactor; >> 146 } >> 147 >> 148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 149 >> 150 G4double G4BraggIonModel::GetParticleCharge(const G4ParticleDefinition* p, >> 151 const G4Material* mat, >> 152 G4double kineticEnergy) >> 153 { >> 154 // this method is called only for ions >> 155 return corr->GetParticleCharge(p,mat,kineticEnergy); >> 156 } >> 157 >> 158 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 159 >> 160 G4double G4BraggIonModel::ComputeCrossSectionPerElectron( >> 161 const G4ParticleDefinition* p, >> 162 G4double kineticEnergy, >> 163 G4double cutEnergy, >> 164 G4double maxKinEnergy) 128 { 165 { 129 // this method is called only for ions, so n << 166 G4double cross = 0.0; 130 if(isAlpha) { return 1.0; } << 167 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy); 131 return G4BraggModel::GetChargeSquareRatio(p, << 168 G4double maxEnergy = std::min(tmax,maxKinEnergy); >> 169 if(cutEnergy < tmax) { >> 170 >> 171 G4double energy = kineticEnergy + mass; >> 172 G4double energy2 = energy*energy; >> 173 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2; >> 174 cross = 1.0/cutEnergy - 1.0/maxEnergy - beta2*log(maxEnergy/cutEnergy)/tmax; >> 175 >> 176 cross *= twopi_mc2_rcl2*chargeSquare/beta2; >> 177 } >> 178 // G4cout << "BR: e= " << kineticEnergy << " tmin= " << cutEnergy >> 179 // << " tmax= " << tmax << " cross= " << cross << G4endl; >> 180 >> 181 return cross; 132 } 182 } 133 183 134 //....oooOO0OOooo........oooOO0OOooo........oo 184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 135 185 136 G4double G4BraggIonModel::ComputeCrossSectionP 186 G4double G4BraggIonModel::ComputeCrossSectionPerAtom( 137 con 187 const G4ParticleDefinition* p, 138 << 188 G4double kineticEnergy, 139 << 189 G4double Z, G4double, 140 190 G4double cutEnergy, 141 191 G4double maxEnergy) 142 { 192 { 143 G4double sigma = << 193 G4double cross = Z*ComputeCrossSectionPerElectron 144 Z*ComputeCrossSectionPerElectron(p,kinEner << 194 (p,kineticEnergy,cutEnergy,maxEnergy); 145 if(isAlpha) { << 195 return cross; 146 sigma *= (HeEffChargeSquare(Z, kinEnergy/C << 147 } << 148 return sigma; << 149 } 196 } 150 197 151 //....oooOO0OOooo........oooOO0OOooo........oo 198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 152 199 153 G4double G4BraggIonModel::CrossSectionPerVolum 200 G4double G4BraggIonModel::CrossSectionPerVolume( 154 con << 201 const G4Material* material, 155 con 202 const G4ParticleDefinition* p, 156 << 203 G4double kineticEnergy, 157 204 G4double cutEnergy, 158 205 G4double maxEnergy) 159 { 206 { 160 G4double sigma = material->GetElectronDensit << 207 G4double eDensity = material->GetElectronDensity(); 161 ComputeCrossSectionPerElectron(p,kinEnergy << 208 G4double cross = eDensity*ComputeCrossSectionPerElectron 162 if(isAlpha) { << 209 (p,kineticEnergy,cutEnergy,maxEnergy); 163 const G4double zeff = material->GetTotNbOf << 210 return cross; 164 material->GetTotNbOfAtomsPerVolume(); << 165 sigma *= (HeEffChargeSquare(zeff, kinEnerg << 166 } << 167 return sigma; << 168 } 211 } 169 212 170 //....oooOO0OOooo........oooOO0OOooo........oo 213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 171 214 172 G4double G4BraggIonModel::ComputeDEDXPerVolume 215 G4double G4BraggIonModel::ComputeDEDXPerVolume(const G4Material* material, 173 << 216 const G4ParticleDefinition* p, 174 << 217 G4double kineticEnergy, 175 << 218 G4double cutEnergy) 176 { << 219 { 177 const G4double tmax = MaxSecondaryEnergy(p, << 220 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy); 178 const G4double tlim = lowestKinEnergy*massRa << 221 G4double tmin = min(cutEnergy, tmax); 179 const G4double tmin = std::max(std::min(cut, << 222 G4double tkin = kineticEnergy/massRate; 180 G4double dedx = 0.0; << 223 G4double dedx = 0.0; >> 224 if(tkin > lowestKinEnergy) dedx = DEDX(material, tkin); >> 225 else dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy); 181 226 182 if(kineticEnergy < tlim) { << 227 if (cutEnergy < tmax) { 183 dedx = HeDEDX(material, tlim)*std::sqrt(ki << 184 } else { << 185 dedx = HeDEDX(material, kineticEnergy); << 186 228 187 if (tmin < tmax) { << 229 G4double tau = kineticEnergy/mass; 188 const G4double tau = kineticEnergy/mass; << 230 G4double gam = tau + 1.0; 189 const G4double x = tmin/tmax; << 231 G4double bg2 = tau * (tau+2.0); 190 << 232 G4double beta2 = bg2/(gam*gam); 191 G4double del = << 233 G4double x = tmin/tmax; 192 (G4Log(x)*(tau + 1.)*(tau + 1.)/(tau * << 234 193 CLHEP::twopi_mc2_rcl2*material->GetElectronD << 235 dedx += (log(x) + (1.0 - x)*beta2) * twopi_mc2_rcl2 194 if(isAlpha) { << 236 * (material->GetElectronDensity())/beta2; 195 const G4double zeff = material->GetTotNbOfEl << 196 material->GetTotNbOfAtomsPerVolume(); << 197 heChargeSquare = HeEffChargeSquare(zeff, kin << 198 del *= heChargeSquare; << 199 } << 200 dedx += del; << 201 } << 202 } 237 } 203 dedx = std::max(dedx, 0.0); << 238 204 /* << 239 // now compute the total ionization loss 205 G4cout << "BraggIon: " << material->GetNam << 240 206 << " E(MeV)=" << kineticEnergy/MeV << 241 if (dedx < 0.0) dedx = 0.0 ; 207 << " Tmin(MeV)=" << tmin << " dedx( << 242 208 << dedx*gram/(MeV*cm2*material->Get << 243 dedx *= chargeSquare; 209 << " q2=" << chargeSquare << G4endl << 244 210 */ << 245 //G4cout << " tkin(MeV) = " << tkin/MeV << " dedx(MeVxcm^2/g) = " >> 246 // << dedx*gram/(MeV*cm2*material->GetDensity()) >> 247 // << " q2 = " << chargeSquare << G4endl; >> 248 211 return dedx; 249 return dedx; 212 } 250 } 213 251 214 //....oooOO0OOooo........oooOO0OOooo........oo 252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 215 253 216 void G4BraggIonModel::CorrectionsAlongStep(con 254 void G4BraggIonModel::CorrectionsAlongStep(const G4MaterialCutsCouple* couple, 217 con << 255 const G4DynamicParticle* dp, 218 con << 256 G4double& eloss, 219 G4d << 257 G4double&, 220 { << 258 G4double length) 221 // no correction for alpha << 259 { 222 if(isAlpha) { return; } << 260 // this method is called only for ions 223 << 224 // no correction at a small step at the last << 225 const G4double preKinEnergy = dp->GetKinetic << 226 if(eloss >= preKinEnergy || eloss < preKinEn << 227 << 228 // corrections only for ions << 229 const G4ParticleDefinition* p = dp->GetDefin 261 const G4ParticleDefinition* p = dp->GetDefinition(); 230 if(p != particle) { SetParticle(p); } << 231 << 232 // effective energy and charge at a step << 233 const G4Material* mat = couple->GetMaterial( 262 const G4Material* mat = couple->GetMaterial(); 234 const G4double e = std::max(preKinEnergy - e << 263 G4double preKinEnergy = dp->GetKineticEnergy(); 235 const G4double q20 = corr->EffectiveChargeSq << 264 G4double e = preKinEnergy - eloss*0.5; 236 const G4double q2 = corr->EffectiveChargeSqu << 265 if(e < 0.0) e = preKinEnergy*0.5; 237 const G4double qfactor = q2/q20; << 266 238 /* << 267 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,e); 239 G4cout << "G4BraggIonModel::CorrectionsAlo << 268 GetModelOfFluctuations()->SetParticleAndCharge(p, q2); 240 << preKinEnergy << " Eeff(MeV)=" << e << 269 eloss *= q2*corr->EffectiveChargeCorrection(p,mat,e)/corrFactor; 241 << " eloss=" << eloss << " elossnew=" << e << 270 242 << " qfactor=" << qfactor << " Qpre=" << q << 271 if(nuclearStopping) { 243 << p->GetParticleName() <<G4endl; << 272 244 */ << 273 G4double nloss = length*corr->NuclearDEDX(p,mat,e,false); 245 eloss *= qfactor; << 274 >> 275 // too big energy loss >> 276 if(eloss + nloss > preKinEnergy) { >> 277 nloss *= (preKinEnergy/(eloss + nloss)); >> 278 eloss = preKinEnergy; >> 279 } else { >> 280 eloss += nloss; >> 281 } >> 282 /* >> 283 G4cout << "G4ionIonisation::CorrectionsAlongStep: e= " << preKinEnergy >> 284 << " de= " << eloss << " NIEL= " << nloss >> 285 << " dynQ= " << dp->GetCharge()/eplus << G4endl; >> 286 */ >> 287 fParticleChange->ProposeNonIonizingEnergyDeposit(nloss); >> 288 } 246 } 289 } 247 290 248 //....oooOO0OOooo........oooOO0OOooo........oo 291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 249 292 250 G4int G4BraggIonModel::HasMaterialForHe(const << 293 void G4BraggIonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp, >> 294 const G4MaterialCutsCouple*, >> 295 const G4DynamicParticle* dp, >> 296 G4double xmin, >> 297 G4double maxEnergy) 251 { 298 { 252 const G4String& chFormula = mat->GetChemical << 299 G4double tmax = MaxSecondaryKinEnergy(dp); 253 if(chFormula.empty()) { return -1; } << 300 G4double xmax = std::min(tmax, maxEnergy); >> 301 if(xmin >= xmax) return; >> 302 >> 303 G4double kineticEnergy = dp->GetKineticEnergy(); >> 304 G4double energy = kineticEnergy + mass; >> 305 G4double energy2 = energy*energy; >> 306 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2; >> 307 G4double grej = 1.0; >> 308 G4double deltaKinEnergy, f; >> 309 >> 310 G4ThreeVector direction = dp->GetMomentumDirection(); >> 311 >> 312 // sampling follows ... >> 313 do { >> 314 G4double q = G4UniformRand(); >> 315 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - q) + xmax*q); >> 316 >> 317 f = 1.0 - beta2*deltaKinEnergy/tmax; >> 318 >> 319 if(f > grej) { >> 320 G4cout << "G4BraggIonModel::SampleSecondary Warning! " >> 321 << "Majorant " << grej << " < " >> 322 << f << " for e= " << deltaKinEnergy >> 323 << G4endl; >> 324 } >> 325 >> 326 } while( grej*G4UniformRand() >= f ); >> 327 >> 328 G4double deltaMomentum = >> 329 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2)); >> 330 G4double totMomentum = energy*sqrt(beta2); >> 331 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) / >> 332 (deltaMomentum * totMomentum); >> 333 if(cost > 1.0) cost = 1.0; >> 334 G4double sint = sqrt((1.0 - cost)*(1.0 + cost)); >> 335 >> 336 G4double phi = twopi * G4UniformRand() ; >> 337 >> 338 G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ; >> 339 deltaDirection.rotateUz(direction); >> 340 >> 341 // create G4DynamicParticle object for delta ray >> 342 G4DynamicParticle* delta = new G4DynamicParticle(theElectron,deltaDirection, >> 343 deltaKinEnergy); >> 344 >> 345 vdp->push_back(delta); >> 346 >> 347 // Change kinematics of primary particle >> 348 kineticEnergy -= deltaKinEnergy; >> 349 G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum; >> 350 finalP = finalP.unit(); >> 351 >> 352 fParticleChange->SetProposedKineticEnergy(kineticEnergy); >> 353 fParticleChange->SetProposedMomentumDirection(finalP); >> 354 } >> 355 >> 356 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 357 >> 358 G4bool G4BraggIonModel::HasMaterial(const G4Material* material) >> 359 { >> 360 const size_t numberOfMolecula = 11 ; >> 361 SetMoleculaNumber(numberOfMolecula) ; >> 362 G4String chFormula = material->GetChemicalFormula() ; 254 363 255 // ICRU Report N49, 1993. Ziegler model for 364 // ICRU Report N49, 1993. Ziegler model for He. 256 << 365 static G4String molName[numberOfMolecula] = { 257 static const G4int numberOfMolecula = 11; << 258 static const G4String molName[numberOfMolecu << 259 "CaF_2", "Cellulose_Nitrate", "LiF", "Po 366 "CaF_2", "Cellulose_Nitrate", "LiF", "Policarbonate", 260 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Po 367 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Polymethly_Methacralate", 261 "Polysterene", "SiO_2", "NaI", "H_2O", 368 "Polysterene", "SiO_2", "NaI", "H_2O", 262 "Graphite" }; << 369 "Graphite" } ; 263 370 264 // Search for the material in the table 371 // Search for the material in the table 265 for (G4int i=0; i<numberOfMolecula; ++i) { << 372 for (size_t i=0; i<numberOfMolecula; i++) { 266 if (chFormula == molName[i]) { << 373 if (chFormula == molName[i]) { 267 return i; << 374 SetMoleculaNumber(i) ; 268 } << 375 return true ; >> 376 } 269 } 377 } 270 return -1; << 378 return false ; 271 } 379 } 272 380 273 //....oooOO0OOooo........oooOO0OOooo........oo 381 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 274 382 275 G4double G4BraggIonModel::HeStoppingPower(cons << 383 G4double G4BraggIonModel::StoppingPower(const G4Material* material, >> 384 G4double kineticEnergy) 276 { 385 { 277 G4double ionloss = 0.0; << 386 G4double ionloss = 0.0 ; 278 if (iMolecula >= 0) { << 387 >> 388 if (iMolecula < 11) { 279 389 280 // The data and the fit from: 390 // The data and the fit from: 281 // ICRU Report N49, 1993. Ziegler's model 391 // ICRU Report N49, 1993. Ziegler's model for alpha 282 // He energy in internal units of parametr 392 // He energy in internal units of parametrisation formula (MeV) 283 // Input scaled energy of a proton or Gene << 284 G4double T = kineticEnergy/(massRate*CLHEP << 285 393 286 static const G4float a[11][5] = { << 394 G4double T = kineticEnergy*rateMassHe2p/MeV ; 287 {9.43672f, 0.54398f, 84.341f, 1.3705f, << 395 288 {67.1503f, 0.41409f, 404.512f, 148.97f, << 396 static G4double a[11][5] = { 289 {5.11203f, 0.453f, 36.718f, 50.6f, << 397 {9.43672, 0.54398, 84.341, 1.3705, 57.422}, 290 {61.793f, 0.48445f, 361.537f, 57.889f, << 398 {67.1503, 0.41409, 404.512, 148.97, 20.99}, 291 {7.83464f, 0.49804f, 160.452f, 3.192f, << 399 {5.11203, 0.453, 36.718, 50.6, 28.058}, 292 {19.729f, 0.52153f, 162.341f, 58.35f, << 400 {61.793, 0.48445, 361.537, 57.889, 50.674}, 293 {26.4648f, 0.50112f, 188.913f, 30.079f, << 401 {7.83464, 0.49804, 160.452, 3.192, 0.71922}, 294 {7.8655f, 0.5205f, 63.96f, 51.32f, << 402 {19.729, 0.52153, 162.341, 58.35, 25.668}, 295 {8.8965f, 0.5148f, 339.36f, 1.7205f, << 403 {26.4648, 0.50112, 188.913, 30.079, 16.509}, 296 {2.959f, 0.53255f, 34.247f, 60.655f, << 404 {7.8655, 0.5205, 63.96, 51.32, 67.775}, 297 {3.80133f, 0.41590f, 12.9966f, 117.83f, << 405 {8.8965, 0.5148, 339.36, 1.7205, 0.70423}, 298 << 406 {2.959, 0.53255, 34.247, 60.655, 15.153}, 299 static const G4double atomicWeight[11] = { << 407 {3.80133, 0.41590, 12.9966, 117.83, 242.28} }; 300 101.96128f, 44.0098f, 16.0426f, 28.0536 << 408 301 104.1512f, 44.665f, 60.0843f, 18.0152 << 409 static G4double atomicWeight[11] = { 302 << 410 101.96128, 44.0098, 16.0426, 28.0536, 42.0804, 303 const G4int i = iMolecula; << 411 104.1512, 44.665, 60.0843, 18.0152, 18.0152, 12.0}; 304 << 412 305 G4double slow = (G4double)(a[i][0]); << 413 G4int i = iMolecula; 306 << 307 G4double x1 = (G4double)(a[i][1]); << 308 G4double x2 = (G4double)(a[i][2]); << 309 G4double x3 = (G4double)(a[i][3]); << 310 G4double x4 = (G4double)(a[i][4]); << 311 414 312 // Free electron gas model 415 // Free electron gas model 313 if ( T < 0.001 ) { 416 if ( T < 0.001 ) { 314 G4double shigh = G4Log( 1.0 + x3*1000.0 << 417 G4double slow = a[i][0] ; >> 418 G4double shigh = log( 1.0 + a[i][3]*1000.0 + a[i][4]*0.001 ) >> 419 * a[i][2]*1000.0 ; 315 ionloss = slow*shigh / (slow + shigh) ; 420 ionloss = slow*shigh / (slow + shigh) ; 316 ionloss *= std::sqrt(T*1000.0) ; << 421 ionloss *= sqrt(T*1000.0) ; 317 422 318 // Main parametrisation 423 // Main parametrisation 319 } else { 424 } else { 320 slow *= G4Exp(G4Log(T*1000.0)*x1) ; << 425 G4double slow = a[i][0] * pow((T*1000.0), a[i][1]) ; 321 G4double shigh = G4Log( 1.0 + x3/T + x4* << 426 G4double shigh = log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ; 322 ionloss = slow*shigh / (slow + shigh) ; 427 ionloss = slow*shigh / (slow + shigh) ; 323 /* 428 /* 324 G4cout << "## " << i << ". T= " << T << 429 G4cout << "## " << i << ". T= " << T << " slow= " << slow 325 << " a0= " << a[i][0] << " a1= " << a << 430 << " a0= " << a[i][0] << " a1= " << a[i][1] 326 << " shigh= " << shigh << 431 << " shigh= " << shigh 327 << " dedx= " << ionloss << " q^2= " < << 432 << " dedx= " << ionloss << " q^2= " << HeEffChargeSquare(z, T*MeV) 328 << G4endl; << 433 << G4endl; 329 */ 434 */ 330 } 435 } 331 ionloss = std::max(ionloss, 0.0) * atomicW << 436 if ( ionloss < 0.0) ionloss = 0.0 ; >> 437 >> 438 // He effective charge >> 439 G4double aa = atomicWeight[iMolecula]; >> 440 ionloss /= (HeEffChargeSquare(0.5*aa, T)*aa); >> 441 >> 442 // pure material (normally not the case for this function) >> 443 } else if(1 == (material->GetNumberOfElements())) { >> 444 G4double z = material->GetZ() ; >> 445 ionloss = ElectronicStoppingPower( z, kineticEnergy ) ; 332 } 446 } >> 447 333 return ionloss; 448 return ionloss; 334 } 449 } 335 450 336 //....oooOO0OOooo........oooOO0OOooo........oo 451 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 337 452 338 G4double G4BraggIonModel::HeElectronicStopping << 453 G4double G4BraggIonModel::ElectronicStoppingPower(G4double z, 339 const G4double kinet << 454 G4double kineticEnergy) const 340 { 455 { 341 G4double ionloss ; 456 G4double ionloss ; 342 G4int i = std::min(z-1, 91); // index of at << 457 G4int i = G4int(z)-1 ; // index of atom 343 //G4cout << "ElectronicStoppingPower z=" << << 458 if(i < 0) i = 0 ; 344 // << " E=" << kineticEnergy << G4endl; << 459 if(i > 91) i = 91 ; >> 460 345 // The data and the fit from: 461 // The data and the fit from: 346 // ICRU Report 49, 1993. Ziegler's type of p 462 // ICRU Report 49, 1993. Ziegler's type of parametrisations. 347 // Proton kinetic energy for parametrisation 463 // Proton kinetic energy for parametrisation (keV/amu) 348 // He energy in internal units of parametris << 464 349 //G4double T = kineticEnergy*rateMassHe2p/CL << 465 // He energy in internal units of parametrisation formula (MeV) 350 G4double T = kineticEnergy/CLHEP::MeV; << 466 G4double T = kineticEnergy*rateMassHe2p/MeV ; 351 << 467 352 static const G4float a[92][5] = { << 468 static G4double a[92][5] = { 353 { 0.35485f, 0.6456f, 6.01525f, 20.8933f, << 469 {0.35485, 0.6456, 6.01525, 20.8933, 4.3515 354 },{ 0.58f, 0.59f, 6.3f, 130.0f, << 470 },{ 0.58, 0.59, 6.3, 130.0, 44.07 355 },{ 1.42f, 0.49f, 12.25f, 32.0f, << 471 },{ 1.42, 0.49, 12.25, 32.0, 9.161 356 },{ 2.206f, 0.51f, 15.32f, 0.25f, << 472 },{ 2.206, 0.51, 15.32, 0.25, 8.995 //Be Ziegler77 357 // },{ 2.1895f, 0.47183,7.2362f, 134 << 473 // },{ 2.1895, 0.47183,7.2362, 134.30, 197.96 //Be from ICRU 358 },{ 3.691f, 0.4128f, 18.48f, 50.72f, << 474 },{ 3.691, 0.4128, 18.48, 50.72, 9.0 359 },{ 3.83523f, 0.42993f,12.6125f, 227.41f, << 475 },{ 3.83523, 0.42993,12.6125, 227.41, 188.97 360 // },{ 1.9259f, 0.5550f, 27.15125f, 26 << 476 },{ 1.9259, 0.5550, 27.15125, 26.0665, 6.2768 361 },{ 1.9259f, 0.5550f, 27.1513f, 26.0665f, << 477 },{ 2.81015, 0.4759, 50.0253, 10.556, 1.0382 362 },{ 2.81015f, 0.4759f, 50.0253f, 10.556f, << 478 },{ 1.533, 0.531, 40.44, 18.41, 2.718 363 },{ 1.533f, 0.531f, 40.44f, 18.41f, << 479 },{ 2.303, 0.4861, 37.01, 37.96, 5.092 364 },{ 2.303f, 0.4861f, 37.01f, 37.96f, << 365 // Z= 11-20 480 // Z= 11-20 366 },{ 9.894f, 0.3081f, 23.65f, 0.384f, << 481 },{ 9.894, 0.3081, 23.65, 0.384, 92.93 367 },{ 4.3f, 0.47f, 34.3f, 3.3f, << 482 },{ 4.3, 0.47, 34.3, 3.3, 12.74 368 },{ 2.5f, 0.625f, 45.7f, 0.1f, << 483 },{ 2.5, 0.625, 45.7, 0.1, 4.359 369 },{ 2.1f, 0.65f, 49.34f, 1.788f, << 484 },{ 2.1, 0.65, 49.34, 1.788, 4.133 370 },{ 1.729f, 0.6562f, 53.41f, 2.405f, << 485 },{ 1.729, 0.6562, 53.41, 2.405, 3.845 371 },{ 1.402f, 0.6791f, 58.98f, 3.528f, << 486 },{ 1.402, 0.6791, 58.98, 3.528, 3.211 372 },{ 1.117f, 0.7044f, 69.69f, 3.705f, << 487 },{ 1.117, 0.7044, 69.69, 3.705, 2.156 373 },{ 2.291f, 0.6284f, 73.88f, 4.478f, << 488 },{ 2.291, 0.6284, 73.88, 4.478, 2.066 374 },{ 8.554f, 0.3817f, 83.61f, 11.84f, << 489 },{ 8.554, 0.3817, 83.61, 11.84, 1.875 375 },{ 6.297f, 0.4622f, 65.39f, 10.14f, << 490 },{ 6.297, 0.4622, 65.39, 10.14, 5.036 376 // Z= 21-30 491 // Z= 21-30 377 },{ 5.307f, 0.4918f, 61.74f, 12.4f, << 492 },{ 5.307, 0.4918, 61.74, 12.4, 6.665 378 },{ 4.71f, 0.5087f, 65.28f, 8.806f, << 493 },{ 4.71, 0.5087, 65.28, 8.806, 5.948 379 },{ 6.151f, 0.4524f, 83.0f, 18.31f, << 494 },{ 6.151, 0.4524, 83.0, 18.31, 2.71 380 },{ 6.57f, 0.4322f, 84.76f, 15.53f, << 495 },{ 6.57, 0.4322, 84.76, 15.53, 2.779 381 },{ 5.738f, 0.4492f, 84.6f, 14.18f, << 496 },{ 5.738, 0.4492, 84.6, 14.18, 3.101 382 },{ 5.013f, 0.4707f, 85.8f, 16.55f, << 497 },{ 5.013, 0.4707, 85.8, 16.55, 3.211 383 },{ 4.32f, 0.4947f, 76.14f, 10.85f, << 498 },{ 4.32, 0.4947, 76.14, 10.85, 5.441 384 },{ 4.652f, 0.4571f, 80.73f, 22.0f, << 499 },{ 4.652, 0.4571, 80.73, 22.0, 4.952 385 },{ 3.114f, 0.5236f, 76.67f, 7.62f, << 500 },{ 3.114, 0.5236, 76.67, 7.62, 6.385 386 },{ 3.114f, 0.5236f, 76.67f, 7.62f, << 501 },{ 3.114, 0.5236, 76.67, 7.62, 7.502 387 // Z= 31-40 502 // Z= 31-40 388 },{ 3.114f, 0.5236f, 76.67f, 7.62f, << 503 },{ 3.114, 0.5236, 76.67, 7.62, 8.514 389 },{ 5.746f, 0.4662f, 79.24f, 1.185f, << 504 },{ 5.746, 0.4662, 79.24, 1.185, 7.993 390 },{ 2.792f, 0.6346f, 106.1f, 0.2986f, << 505 },{ 2.792, 0.6346, 106.1, 0.2986, 2.331 391 },{ 4.667f, 0.5095f, 124.3f, 2.102f, << 506 },{ 4.667, 0.5095, 124.3, 2.102, 1.667 392 },{ 2.44f, 0.6346f, 105.0f, 0.83f, << 507 },{ 2.44, 0.6346, 105.0, 0.83, 2.851 393 },{ 1.413f, 0.7377f, 147.9f, 1.466f, << 508 },{ 1.413, 0.7377, 147.9, 1.466, 1.016 394 },{ 11.72f, 0.3826f, 102.8f, 9.231f, << 509 },{ 11.72, 0.3826, 102.8, 9.231, 4.371 395 },{ 7.126f, 0.4804f, 119.3f, 5.784f, << 510 },{ 7.126, 0.4804, 119.3, 5.784, 2.454 396 },{ 11.61f, 0.3955f, 146.7f, 7.031f, << 511 },{ 11.61, 0.3955, 146.7, 7.031, 1.423 397 },{ 10.99f, 0.41f, 163.9f, 7.1f, << 512 },{ 10.99, 0.41, 163.9, 7.1, 1.052 398 // Z= 41-50 513 // Z= 41-50 399 },{ 9.241f, 0.4275f, 163.1f, 7.954f, << 514 },{ 9.241, 0.4275, 163.1, 7.954, 1.102 400 },{ 9.276f, 0.418f, 157.1f, 8.038f, << 515 },{ 9.276, 0.418, 157.1, 8.038, 1.29 401 },{ 3.999f, 0.6152f, 97.6f, 1.297f, << 516 },{ 3.999, 0.6152, 97.6, 1.297, 5.792 402 },{ 4.306f, 0.5658f, 97.99f, 5.514f, << 517 },{ 4.306, 0.5658, 97.99, 5.514, 5.754 403 },{ 3.615f, 0.6197f, 86.26f, 0.333f, << 518 },{ 3.615, 0.6197, 86.26, 0.333, 8.689 404 },{ 5.8f, 0.49f, 147.2f, 6.903f, << 519 },{ 5.8, 0.49, 147.2, 6.903, 1.289 405 },{ 5.6f, 0.49f, 130.0f, 10.0f, << 520 },{ 5.6, 0.49, 130.0, 10.0, 2.844 406 },{ 3.55f, 0.6068f, 124.7f, 1.112f, << 521 },{ 3.55, 0.6068, 124.7, 1.112, 3.119 407 },{ 3.6f, 0.62f, 105.8f, 0.1692f, << 522 },{ 3.6, 0.62, 105.8, 0.1692, 6.026 408 },{ 5.4f, 0.53f, 103.1f, 3.931f, << 523 },{ 5.4, 0.53, 103.1, 3.931, 7.767 409 // Z= 51-60 524 // Z= 51-60 410 },{ 3.97f, 0.6459f, 131.8f, 0.2233f, << 525 },{ 3.97, 0.6459, 131.8, 0.2233, 2.723 411 },{ 3.65f, 0.64f, 126.8f, 0.6834f, << 526 },{ 3.65, 0.64, 126.8, 0.6834, 3.411 412 },{ 3.118f, 0.6519f, 164.9f, 1.208f, << 527 },{ 3.118, 0.6519, 164.9, 1.208, 1.51 413 },{ 3.949f, 0.6209f, 200.5f, 1.878f, << 528 },{ 3.949, 0.6209, 200.5, 1.878, 0.9126 414 },{ 14.4f, 0.3923f, 152.5f, 8.354f, << 529 },{ 14.4, 0.3923, 152.5, 8.354, 2.597 415 },{ 10.99f, 0.4599f, 138.4f, 4.811f, << 530 },{ 10.99, 0.4599, 138.4, 4.811, 3.726 416 },{ 16.6f, 0.3773f, 224.1f, 6.28f, << 531 },{ 16.6, 0.3773, 224.1, 6.28, 0.9121 417 },{ 10.54f, 0.4533f, 159.3f, 4.832f, << 532 },{ 10.54, 0.4533, 159.3, 4.832, 2.529 418 },{ 10.33f, 0.4502f, 162.0f, 5.132f, << 533 },{ 10.33, 0.4502, 162.0, 5.132, 2.444 419 },{ 10.15f, 0.4471f, 165.6f, 5.378f, << 534 },{ 10.15, 0.4471, 165.6, 5.378, 2.328 420 // Z= 61-70 535 // Z= 61-70 421 },{ 9.976f, 0.4439f, 168.0f, 5.721f, << 536 },{ 9.976, 0.4439, 168.0, 5.721, 2.258 422 },{ 9.804f, 0.4408f, 176.2f, 5.675f, << 537 },{ 9.804, 0.4408, 176.2, 5.675, 1.997 423 },{ 14.22f, 0.363f, 228.4f, 7.024f, << 538 },{ 14.22, 0.363, 228.4, 7.024, 1.016 424 },{ 9.952f, 0.4318f, 233.5f, 5.065f, << 539 },{ 9.952, 0.4318, 233.5, 5.065, 0.9244 425 },{ 9.272f, 0.4345f, 210.0f, 4.911f, << 540 },{ 9.272, 0.4345, 210.0, 4.911, 1.258 426 },{ 10.13f, 0.4146f, 225.7f, 5.525f, << 541 },{ 10.13, 0.4146, 225.7, 5.525, 1.055 427 },{ 8.949f, 0.4304f, 213.3f, 5.071f, << 542 },{ 8.949, 0.4304, 213.3, 5.071, 1.221 428 },{ 11.94f, 0.3783f, 247.2f, 6.655f, << 543 },{ 11.94, 0.3783, 247.2, 6.655, 0.849 429 },{ 8.472f, 0.4405f, 195.5f, 4.051f, << 544 },{ 8.472, 0.4405, 195.5, 4.051, 1.604 430 },{ 8.301f, 0.4399f, 203.7f, 3.667f, << 545 },{ 8.301, 0.4399, 203.7, 3.667, 1.459 431 // Z= 71-80 546 // Z= 71-80 432 },{ 6.567f, 0.4858f, 193.0f, 2.65f, << 547 },{ 6.567, 0.4858, 193.0, 2.65, 1.66 433 },{ 5.951f, 0.5016f, 196.1f, 2.662f, << 548 },{ 5.951, 0.5016, 196.1, 2.662, 1.589 434 },{ 7.495f, 0.4523f, 251.4f, 3.433f, << 549 },{ 7.495, 0.4523, 251.4, 3.433, 0.8619 435 },{ 6.335f, 0.4825f, 255.1f, 2.834f, << 550 },{ 6.335, 0.4825, 255.1, 2.834, 0.8228 436 },{ 4.314f, 0.5558f, 214.8f, 2.354f, << 551 },{ 4.314, 0.5558, 214.8, 2.354, 1.263 437 },{ 4.02f, 0.5681f, 219.9f, 2.402f, << 552 },{ 4.02, 0.5681, 219.9, 2.402, 1.191 438 },{ 3.836f, 0.5765f, 210.2f, 2.742f, << 553 },{ 3.836, 0.5765, 210.2, 2.742, 1.305 439 },{ 4.68f, 0.5247f, 244.7f, 2.749f, << 554 },{ 4.68, 0.5247, 244.7, 2.749, 0.8962 440 },{ 2.892f, 0.6204f, 208.6f, 2.415f, << 555 },{ 2.892, 0.6204, 208.6, 2.415, 1.416 //Au Z77 441 // },{ 3.223f, 0.5883f, 232.7f, 2.9 << 556 // },{ 3.223, 0.5883, 232.7, 2.954, 1.05 //Au ICRU 442 },{ 2.892f, 0.6204f, 208.6f, 2.415f, << 557 },{ 2.892, 0.6204, 208.6, 2.415, 1.416 443 // Z= 81-90 558 // Z= 81-90 444 },{ 4.728f, 0.5522f, 217.0f, 3.091f, << 559 },{ 4.728, 0.5522, 217.0, 3.091, 1.386 445 },{ 6.18f, 0.52f, 170.0f, 4.0f, << 560 },{ 6.18, 0.52, 170.0, 4.0, 3.224 446 },{ 9.0f, 0.47f, 198.0f, 3.8f, << 561 },{ 9.0, 0.47, 198.0, 3.8, 2.032 447 },{ 2.324f, 0.6997f, 216.0f, 1.599f, << 562 },{ 2.324, 0.6997, 216.0, 1.599, 1.399 448 },{ 1.961f, 0.7286f, 223.0f, 1.621f, << 563 },{ 1.961, 0.7286, 223.0, 1.621, 1.296 449 },{ 1.75f, 0.7427f, 350.1f, 0.9789f, << 564 },{ 1.75, 0.7427, 350.1, 0.9789, 0.5507 450 },{ 10.31f, 0.4613f, 261.2f, 4.738f, << 565 },{ 10.31, 0.4613, 261.2, 4.738, 0.9899 451 },{ 7.962f, 0.519f, 235.7f, 4.347f, << 566 },{ 7.962, 0.519, 235.7, 4.347, 1.313 452 },{ 6.227f, 0.5645f, 231.9f, 3.961f, << 567 },{ 6.227, 0.5645, 231.9, 3.961, 1.379 453 },{ 5.246f, 0.5947f, 228.6f, 4.027f, << 568 },{ 5.246, 0.5947, 228.6, 4.027, 1.432 454 // Z= 91-92 569 // Z= 91-92 455 },{ 5.408f, 0.5811f, 235.7f, 3.961f, << 570 },{ 5.408, 0.5811, 235.7, 3.961, 1.358 456 },{ 5.218f, 0.5828f, 245.0f, 3.838f, << 571 },{ 5.218, 0.5828, 245.0, 3.838, 1.25} 457 }; 572 }; 458 573 459 G4double slow = (G4double)(a[i][0]); << 460 << 461 G4double x1 = (G4double)(a[i][1]); << 462 G4double x2 = (G4double)(a[i][2]); << 463 G4double x3 = (G4double)(a[i][3]); << 464 G4double x4 = (G4double)(a[i][4]); << 465 << 466 // Free electron gas model 574 // Free electron gas model 467 if ( T < 0.001 ) { 575 if ( T < 0.001 ) { 468 G4double shigh = G4Log( 1.0 + x3*1000.0 + << 576 G4double slow = a[i][0] ; 469 ionloss = slow*shigh*std::sqrt(T*1000.0) << 577 G4double shigh = log( 1.0 + a[i][3]*1000.0 + a[i][4]*0.001 ) >> 578 * a[i][2]*1000.0 ; >> 579 ionloss = slow*shigh / (slow + shigh) ; >> 580 ionloss *= sqrt(T*1000.0) ; 470 581 471 // Main parametrisation 582 // Main parametrisation 472 } else { 583 } else { 473 slow *= G4Exp(G4Log(T*1000.0)*x1); << 584 G4double slow = a[i][0] * pow((T*1000.0), a[i][1]) ; 474 G4double shigh = G4Log( 1.0 + x3/T + x4*T << 585 G4double shigh = log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ; 475 ionloss = slow*shigh / (slow + shigh) ; 586 ionloss = slow*shigh / (slow + shigh) ; 476 /* 587 /* 477 G4cout << "## " << i << ". T= " << T << " 588 G4cout << "## " << i << ". T= " << T << " slow= " << slow 478 << " a0= " << a[i][0] << " a1= " << 589 << " a0= " << a[i][0] << " a1= " << a[i][1] 479 << " shigh= " << shigh 590 << " shigh= " << shigh 480 << " dedx= " << ionloss << " q^2= " << 591 << " dedx= " << ionloss << " q^2= " << HeEffChargeSquare(z, T*MeV) 481 << G4endl; << 592 << G4endl; 482 */ 593 */ 483 } 594 } 484 ionloss = std::max(ionloss, 0.0); << 595 if ( ionloss < 0.0) ionloss = 0.0 ; >> 596 >> 597 // He effective charge >> 598 ionloss /= HeEffChargeSquare(z, T); >> 599 485 return ionloss; 600 return ionloss; 486 } 601 } 487 602 488 //....oooOO0OOooo........oooOO0OOooo........oo 603 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 489 604 490 G4double G4BraggIonModel::HeDEDX(const G4Mater << 605 G4double G4BraggIonModel::DEDX(const G4Material* material, 491 const G4double << 606 G4double kineticEnergy) 492 { 607 { 493 // aEnergy is energy of alpha << 494 G4double eloss = 0.0; 608 G4double eloss = 0.0; 495 // check DB << 609 const G4int numberOfElements = material->GetNumberOfElements(); 496 if(material != currentMaterial) { << 497 currentMaterial = material; << 498 baseMaterial = material->GetBaseMaterial() << 499 ? material->GetBaseMaterial() : material << 500 iPSTAR = -1; << 501 iASTAR = -1; << 502 iMolecula = -1; << 503 iICRU90 = (nullptr != fICRU90) ? fICRU90-> << 504 << 505 if(iICRU90 < 0) { << 506 if(isAlpha) { << 507 iASTAR = fASTAR->GetIndex(baseMaterial); << 508 if(iASTAR < 0) { iMolecula = HasMaterialForH << 509 } else { << 510 iPSTAR = fPSTAR->GetIndex(baseMaterial); << 511 } << 512 } << 513 /* << 514 G4cout << "%%% " <<material->GetName() << << 515 << iMolecula << " iASTAR= " << iAS << 516 << " iICRU90= " << iICRU90<< G4end << 517 */ << 518 } << 519 // ICRU90 << 520 if(iICRU90 >= 0) { << 521 eloss = (isAlpha) << 522 ? fICRU90->GetElectronicDEDXforAlpha(iIC << 523 : fICRU90->GetElectronicDEDXforProton(iI << 524 if(eloss > 0.0) { return eloss*material->G << 525 } << 526 // PSTAR parameterisation << 527 if( iPSTAR >= 0 ) { << 528 return fPSTAR->GetElectronicDEDX(iPSTAR, a << 529 *material->GetDensity(); << 530 } << 531 // ASTAR << 532 if( iASTAR >= 0 ) { << 533 eloss = fASTAR->GetElectronicDEDX(iASTAR, << 534 /* << 535 G4cout << "ASTAR: E=" << aEnergy << 536 << " dedx=" << eloss*material->GetDensity << 537 << " " << particle->GetParticleName() << << 538 */ << 539 if(eloss > 0.0) { return eloss*material->G << 540 } << 541 << 542 const std::size_t numberOfElements = materia << 543 const G4ElementVector* theElmVector = materi << 544 const G4double* theAtomicNumDensityVector = 610 const G4double* theAtomicNumDensityVector = 545 material->GetAtomicNumDensityVector(); << 611 material->GetAtomicNumDensityVector(); >> 612 >> 613 // compaund material with parametrisation >> 614 G4int iNist = astar.GetIndex(material); 546 615 547 // molecular data use proton stopping power << 616 if( iNist >= 0 ) { 548 // element data from ICRU49 include data for << 617 G4double T = kineticEnergy*rateMassHe2p; 549 if(iMolecula >= 0) { << 618 return astar.GetElectronicDEDX(iNist, T)*material->GetDensity()/ 550 const G4double zeff = material->GetTotNbOf << 619 HeEffChargeSquare(astar.GetEffectiveZ(iNist), T/MeV); 551 material->GetTotNbOfAtomsPerVolume(); << 552 heChargeSquare = HeEffChargeSquare(zeff, a << 553 eloss = HeStoppingPower(aEnergy)*heChargeS << 554 620 555 // pure material << 621 } else if( HasMaterial(material) ) { >> 622 >> 623 eloss = StoppingPower(material, kineticEnergy)* >> 624 material->GetDensity()/amu; >> 625 >> 626 // pure material 556 } else if(1 == numberOfElements) { 627 } else if(1 == numberOfElements) { 557 628 558 const G4Element* element = (*theElmVector) << 629 G4double z = material->GetZ(); 559 eloss = HeElectronicStoppingPower(element- << 630 eloss = ElectronicStoppingPower(z, kineticEnergy) 560 * (material->GetTotNbOfAtomsPerVolume()) << 631 * (material->GetTotNbOfAtomsPerVolume()); 561 632 562 // Brugg's rule calculation 633 // Brugg's rule calculation 563 } else { 634 } else { >> 635 const G4ElementVector* theElementVector = >> 636 material->GetElementVector() ; >> 637 564 // loop for the elements in the material 638 // loop for the elements in the material 565 for (std::size_t i=0; i<numberOfElements; << 639 for (G4int i=0; i<numberOfElements; i++) 566 const G4Element* element = (*theElmVecto << 640 { 567 eloss += HeElectronicStoppingPower(eleme << 641 const G4Element* element = (*theElementVector)[i] ; 568 * theAtomicNumDensityVector[i]; << 642 eloss += ElectronicStoppingPower(element->GetZ(), kineticEnergy) >> 643 * theAtomicNumDensityVector[i]; 569 } 644 } 570 } 645 } 571 return eloss*theZieglerFactor; 646 return eloss*theZieglerFactor; 572 } 647 } 573 648 574 //....oooOO0OOooo........oooOO0OOooo........oo 649 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 575 650 576 G4double << 651 G4double G4BraggIonModel::HeEffChargeSquare(G4double z, 577 G4BraggIonModel::HeEffChargeSquare(const G4dou << 652 G4double kinEnergyHeInMeV) const 578 const G4dou << 579 { 653 { 580 // The aproximation of He effective charge f 654 // The aproximation of He effective charge from: 581 // J.F.Ziegler, J.P. Biersack, U. Littmark 655 // J.F.Ziegler, J.P. Biersack, U. Littmark 582 // The Stopping and Range of Ions in Matter, 656 // The Stopping and Range of Ions in Matter, 583 // Vol.1, Pergamon Press, 1985 657 // Vol.1, Pergamon Press, 1985 584 658 585 static const G4double c[6] = {0.2865, 0.126 << 659 static G4double c[6] = {0.2865, 0.1266, -0.001429, 586 0.02402,-0.011 << 660 0.02402,-0.01135, 0.001475}; 587 661 588 G4double e = std::max(0.0, G4Log(kinEnergyHe << 662 G4double e = std::max(0.0,std::log(kinEnergyHeInMeV*massFactor)); 589 G4double x = c[0] ; 663 G4double x = c[0] ; 590 G4double y = 1.0 ; 664 G4double y = 1.0 ; 591 for (G4int i=1; i<6; ++i) { << 665 for (G4int i=1; i<6; i++) { 592 y *= e; << 666 y *= e ; 593 x += y * c[i]; << 667 x += y * c[i] ; 594 } 668 } 595 669 596 G4double w = 7.6 - e ; 670 G4double w = 7.6 - e ; 597 w = 1.0 + (0.007 + 0.00005*z) * G4Exp( -w*w << 671 w = 1.0 + (0.007 + 0.00005*z) * exp( -w*w ) ; 598 w = 4.0 * (1.0 - G4Exp(-x)) * w * w ; << 672 w = 4.0 * (1.0 - exp(-x)) * w * w ; 599 673 600 return w; 674 return w; 601 } 675 } 602 676 603 //....oooOO0OOooo........oooOO0OOooo........oo 677 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 604 678 605 679