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