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: G4BetheBlochModel.cc,v 1.24 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 header file 31 // GEANT4 Class header file 29 // 32 // 30 // 33 // 31 // File name: G4BetheBlochModel 34 // File name: G4BetheBlochModel 32 // 35 // 33 // Author: Vladimir Ivanchenko on base 36 // Author: Vladimir Ivanchenko on base of Laszlo Urban code 34 // 37 // 35 // Creation date: 03.01.2002 38 // Creation date: 03.01.2002 36 // 39 // 37 // Modifications: 40 // Modifications: 38 // 41 // 39 // 04-12-02 Fix problem of G4DynamicParticle c 42 // 04-12-02 Fix problem of G4DynamicParticle constructor (V.Ivanchenko) 40 // 23-12-02 Change interface in order to move 43 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko) 41 // 27-01-03 Make models region aware (V.Ivanch 44 // 27-01-03 Make models region aware (V.Ivanchenko) 42 // 13-02-03 Add name (V.Ivanchenko) 45 // 13-02-03 Add name (V.Ivanchenko) 43 // 24-03-05 Add G4EmCorrections (V.Ivanchenko) 46 // 24-03-05 Add G4EmCorrections (V.Ivanchenko) 44 // 11-04-05 Major optimisation of internal int 47 // 11-04-05 Major optimisation of internal interfaces (V.Ivanchenko) 45 // 11-02-06 ComputeCrossSectionPerElectron, Co 48 // 11-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma) 46 // 12-02-06 move G4LossTableManager::Instance( 49 // 12-02-06 move G4LossTableManager::Instance()->EmCorrections() 47 // in constructor (mma) 50 // in constructor (mma) 48 // 12-08-08 Added methods GetParticleCharge, G 51 // 12-08-08 Added methods GetParticleCharge, GetChargeSquareRatio, 49 // CorrectionsAlongStep needed for io 52 // CorrectionsAlongStep needed for ions(V.Ivanchenko) 50 // 53 // 51 // ------------------------------------------- 54 // ------------------------------------------------------------------- 52 // 55 // 53 56 >> 57 54 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 55 //....oooOO0OOooo........oooOO0OOooo........oo 59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 56 60 57 #include "G4BetheBlochModel.hh" 61 #include "G4BetheBlochModel.hh" 58 #include "Randomize.hh" 62 #include "Randomize.hh" 59 #include "G4PhysicalConstants.hh" << 60 #include "G4SystemOfUnits.hh" << 61 #include "G4NistManager.hh" << 62 #include "G4Electron.hh" 63 #include "G4Electron.hh" 63 #include "G4LossTableManager.hh" 64 #include "G4LossTableManager.hh" 64 #include "G4EmCorrections.hh" 65 #include "G4EmCorrections.hh" 65 #include "G4EmParameters.hh" << 66 #include "G4ParticleChangeForLoss.hh" 66 #include "G4ParticleChangeForLoss.hh" 67 #include "G4ICRU90StoppingData.hh" << 67 #include "G4NistManager.hh" 68 #include "G4Log.hh" << 69 #include "G4DeltaAngle.hh" << 70 #include <vector> << 71 68 72 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 73 70 74 G4BetheBlochModel::G4BetheBlochModel(const G4P << 71 using namespace std; >> 72 >> 73 G4BetheBlochModel::G4BetheBlochModel(const G4ParticleDefinition* p, 75 const G4S 74 const G4String& nam) 76 : G4VEmModel(nam), 75 : G4VEmModel(nam), 77 twoln10(2.0*G4Log(10.0)), << 76 particle(0), 78 fAlphaTlimit(1*CLHEP::GeV), << 77 tlimit(DBL_MAX), 79 fProtonTlimit(10*CLHEP::GeV) << 78 twoln10(2.0*log(10.0)), >> 79 bg2lim(0.0169), >> 80 taulim(8.4146e-3), >> 81 isIon(false), >> 82 isInitialised(false) 80 { 83 { >> 84 fParticleChange = 0; >> 85 if(p) SetParticle(p); 81 theElectron = G4Electron::Electron(); 86 theElectron = G4Electron::Electron(); 82 corr = G4LossTableManager::Instance()->EmCor 87 corr = G4LossTableManager::Instance()->EmCorrections(); 83 nist = G4NistManager::Instance(); 88 nist = G4NistManager::Instance(); 84 SetLowEnergyLimit(2.0*CLHEP::MeV); << 89 SetLowEnergyLimit(2.0*MeV); 85 } 90 } 86 91 87 //....oooOO0OOooo........oooOO0OOooo........oo 92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 88 93 89 G4BetheBlochModel::~G4BetheBlochModel() = defa << 94 G4BetheBlochModel::~G4BetheBlochModel() >> 95 {} >> 96 >> 97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 98 >> 99 G4double G4BetheBlochModel::MinEnergyCut(const G4ParticleDefinition*, >> 100 const G4MaterialCutsCouple* couple) >> 101 { >> 102 return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy(); >> 103 } 90 104 91 //....oooOO0OOooo........oooOO0OOooo........oo 105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 92 106 93 void G4BetheBlochModel::Initialise(const G4Par 107 void G4BetheBlochModel::Initialise(const G4ParticleDefinition* p, 94 const G4Dat 108 const G4DataVector&) 95 { 109 { 96 if(p != particle) { SetupParameters(p); } << 110 if (!particle) SetParticle(p); 97 111 98 // always false before the run << 112 corrFactor = chargeSquare; 99 SetDeexcitationFlag(false); << 100 113 101 // initialisation once << 114 if(!isInitialised) { 102 if(nullptr == fParticleChange) { << 115 isInitialised = true; 103 const G4String& pname = particle->GetParti << 104 if(G4EmParameters::Instance()->UseICRU90Da << 105 (pname == "proton" || pname == "Generic << 106 fICRU90 = nist->GetICRU90StoppingData(); << 107 } << 108 if (pname == "GenericIon") { << 109 isIon = true; << 110 } else if (pname == "alpha") { << 111 isAlpha = true; << 112 } else if (particle->GetPDGCharge() > 1.1* << 113 isIon = true; << 114 } << 115 116 116 fParticleChange = GetParticleChangeForLoss << 117 if(!fParticleChange) { 117 if(UseAngularGeneratorFlag() && nullptr == << 118 if (pParticleChange) { 118 SetAngularDistribution(new G4DeltaAngle( << 119 fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*> >> 120 (pParticleChange); >> 121 } else { >> 122 fParticleChange = new G4ParticleChangeForLoss(); >> 123 } 119 } 124 } 120 } 125 } 121 // initialisation for each new run << 122 if(IsMaster() && nullptr != fICRU90) { << 123 fICRU90->Initialise(); << 124 } << 125 } 126 } 126 127 127 //....oooOO0OOooo........oooOO0OOooo........oo << 128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 128 129 129 G4double G4BetheBlochModel::GetChargeSquareRat << 130 void G4BetheBlochModel::SetParticle(const G4ParticleDefinition* p) 130 << 131 << 132 { 131 { 133 // this method is called only for ions, so n << 132 if(particle != p) { 134 if(isAlpha) { return 1.0; } << 133 particle = p; 135 chargeSquare = corr->EffectiveChargeSquareRa << 134 G4String pname = particle->GetParticleName(); 136 return chargeSquare; << 135 if (particle->GetParticleType() == "nucleus" && >> 136 pname != "deuteron" && pname != "triton") { >> 137 isIon = true; >> 138 } >> 139 >> 140 mass = particle->GetPDGMass(); >> 141 spin = particle->GetPDGSpin(); >> 142 G4double q = particle->GetPDGCharge()/eplus; >> 143 chargeSquare = q*q; >> 144 ratio = electron_mass_c2/mass; >> 145 G4double magmom = particle->GetPDGMagneticMoment() >> 146 *mass/(0.5*eplus*hbar_Planck*c_squared); >> 147 magMoment2 = magmom*magmom - 1.0; >> 148 formfact = 0.0; >> 149 if(particle->GetLeptonNumber() == 0) { >> 150 G4double x = 0.8426*GeV; >> 151 if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;} >> 152 else if(mass > GeV) { >> 153 x /= nist->GetZ13(mass/proton_mass_c2); >> 154 // tlimit = 51.2*GeV*A13[iz]*A13[iz]; >> 155 } >> 156 formfact = 2.0*electron_mass_c2/(x*x); >> 157 tlimit = 2.0/formfact; >> 158 } >> 159 } 137 } 160 } 138 161 139 //....oooOO0OOooo........oooOO0OOooo........oo 162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 140 163 141 G4double G4BetheBlochModel::GetParticleCharge( << 164 G4double G4BetheBlochModel::GetChargeSquareRatio(const G4ParticleDefinition* p, 142 << 165 const G4Material* mat, 143 << 166 G4double kineticEnergy) 144 { << 145 // this method is called only for ions, so n << 146 return corr->GetParticleCharge(p, mat, kinet << 147 } << 148 << 149 //....oooOO0OOooo........oooOO0OOooo........oo << 150 << 151 void G4BetheBlochModel::SetupParameters(const << 152 { 167 { 153 particle = p; << 168 // this method is called only for ions 154 mass = particle->GetPDGMass(); << 169 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy); 155 spin = particle->GetPDGSpin(); << 170 corrFactor = q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy); 156 G4double q = particle->GetPDGCharge()*invepl << 171 return corrFactor; 157 chargeSquare = q*q; << 158 ratio = electron_mass_c2/mass; << 159 constexpr G4double aMag = 1./(0.5*eplus*CLHE << 160 G4double magmom = particle->GetPDGMagneticMo << 161 magMoment2 = magmom*magmom - 1.0; << 162 formfact = 0.0; << 163 tlimit = DBL_MAX; << 164 if(particle->GetLeptonNumber() == 0) { << 165 G4double x = 0.8426*CLHEP::GeV; << 166 if(spin == 0.0 && mass < CLHEP::GeV) { x = << 167 else if (mass > CLHEP::GeV) { << 168 G4int iz = G4lrint(std::abs(q)); << 169 if(iz > 1) { x /= nist->GetA27(iz); } << 170 } << 171 formfact = 2.0*CLHEP::electron_mass_c2/(x* << 172 tlimit = 2.0/formfact; << 173 } << 174 } 172 } 175 173 176 //....oooOO0OOooo........oooOO0OOooo........oo 174 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 177 175 178 G4double G4BetheBlochModel::MinEnergyCut(const << 176 G4double G4BetheBlochModel::GetParticleCharge(const G4ParticleDefinition* p, 179 const << 177 const G4Material* mat, >> 178 G4double kineticEnergy) 180 { 179 { 181 return couple->GetMaterial()->GetIonisation( << 180 // this method is called only for ions >> 181 return corr->GetParticleCharge(p,mat,kineticEnergy); 182 } 182 } 183 183 184 //....oooOO0OOooo........oooOO0OOooo........oo 184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 185 185 186 G4double 186 G4double 187 G4BetheBlochModel::ComputeCrossSectionPerElect 187 G4BetheBlochModel::ComputeCrossSectionPerElectron(const G4ParticleDefinition* p, 188 << 188 G4double kineticEnergy, 189 << 189 G4double cutEnergy, 190 << 190 G4double maxKinEnergy) 191 { 191 { 192 G4double cross = 0.0; 192 G4double cross = 0.0; 193 const G4double tmax = MaxSecondaryEnergy(p, << 193 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy); 194 const G4double cutEnergy = std::min(std::min << 194 G4double maxEnergy = min(tmax,maxKinEnergy); 195 const G4double maxEnergy = std::min(tmax, ma << 196 if(cutEnergy < maxEnergy) { 195 if(cutEnergy < maxEnergy) { 197 196 198 G4double totEnergy = kineticEnergy + mass; 197 G4double totEnergy = kineticEnergy + mass; 199 G4double energy2 = totEnergy*totEnergy; 198 G4double energy2 = totEnergy*totEnergy; 200 G4double beta2 = kineticEnergy*(kineti 199 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2; 201 200 202 cross = (maxEnergy - cutEnergy)/(cutEnergy << 201 cross = 1.0/cutEnergy - 1.0/maxEnergy 203 - beta2*G4Log(maxEnergy/cutEnergy)/tmax; << 202 - beta2*log(maxEnergy/cutEnergy)/tmax; 204 203 205 // +term for spin=1/2 particle 204 // +term for spin=1/2 particle 206 if( 0.0 < spin ) { cross += 0.5*(maxEnergy << 205 if( 0.5 == spin ) cross += 0.5*(maxEnergy - cutEnergy)/energy2; >> 206 >> 207 // High order correction different for hadrons and ions >> 208 // nevetheless they are applied to reduce high energy transfers >> 209 // if(!isIon) >> 210 //cross += corr->FiniteSizeCorrectionXS(p,currentMaterial, >> 211 // kineticEnergy,cutEnergy); 207 212 208 cross *= CLHEP::twopi_mc2_rcl2*chargeSquar << 213 cross *= twopi_mc2_rcl2*chargeSquare/beta2; 209 } 214 } 210 215 211 // G4cout << "BB: e= " << kineticEnergy << 216 // G4cout << "BB: e= " << kineticEnergy << " tmin= " << cutEnergy 212 // << " tmax= " << tmax << " cross= 217 // << " tmax= " << tmax << " cross= " << cross << G4endl; 213 218 214 return cross; 219 return cross; 215 } 220 } 216 221 217 //....oooOO0OOooo........oooOO0OOooo........oo 222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 218 223 219 G4double G4BetheBlochModel::ComputeCrossSectio 224 G4double G4BetheBlochModel::ComputeCrossSectionPerAtom( 220 con 225 const G4ParticleDefinition* p, 221 << 226 G4double kineticEnergy, 222 << 227 G4double Z, G4double, 223 228 G4double cutEnergy, 224 229 G4double maxEnergy) 225 { 230 { 226 return Z*ComputeCrossSectionPerElectron(p,ki << 231 G4double cross = Z*ComputeCrossSectionPerElectron >> 232 (p,kineticEnergy,cutEnergy,maxEnergy); >> 233 return cross; 227 } 234 } 228 235 229 //....oooOO0OOooo........oooOO0OOooo........oo 236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 230 237 231 G4double G4BetheBlochModel::CrossSectionPerVol 238 G4double G4BetheBlochModel::CrossSectionPerVolume( 232 con << 239 const G4Material* material, 233 con 240 const G4ParticleDefinition* p, 234 << 241 G4double kineticEnergy, 235 242 G4double cutEnergy, 236 243 G4double maxEnergy) 237 { 244 { 238 G4double sigma = mat->GetElectronDensity() << 245 currentMaterial = material; 239 *ComputeCrossSectionPerElectron(p,kinEnerg << 246 G4double eDensity = material->GetElectronDensity(); 240 if(isAlpha) { << 247 G4double cross = eDensity*ComputeCrossSectionPerElectron 241 sigma *= corr->EffectiveChargeSquareRatio( << 248 (p,kineticEnergy,cutEnergy,maxEnergy); 242 } << 249 return cross; 243 return sigma; << 244 } 250 } 245 251 246 //....oooOO0OOooo........oooOO0OOooo........oo 252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 247 253 248 G4double G4BetheBlochModel::ComputeDEDXPerVolu 254 G4double G4BetheBlochModel::ComputeDEDXPerVolume(const G4Material* material, 249 << 255 const G4ParticleDefinition* p, 250 << 256 G4double kineticEnergy, 251 << 257 G4double cut) 252 { 258 { 253 const G4double tmax = MaxSecondaryEnergy(p, << 259 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy); 254 // projectile formfactor limit energy loss << 260 G4double cutEnergy = min(cut,tmax); 255 const G4double cutEnergy = std::min(std::min << 256 261 257 G4double tau = kineticEnergy/mass; 262 G4double tau = kineticEnergy/mass; 258 G4double gam = tau + 1.0; 263 G4double gam = tau + 1.0; 259 G4double bg2 = tau * (tau+2.0); 264 G4double bg2 = tau * (tau+2.0); 260 G4double beta2 = bg2/(gam*gam); 265 G4double beta2 = bg2/(gam*gam); 261 G4double xc = cutEnergy/tmax; << 262 266 263 G4double eexc = material->GetIonisation()-> 267 G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy(); 264 G4double eexc2 = eexc*eexc; 268 G4double eexc2 = eexc*eexc; >> 269 G4double cden = material->GetIonisation()->GetCdensity(); >> 270 G4double mden = material->GetIonisation()->GetMdensity(); >> 271 G4double aden = material->GetIonisation()->GetAdensity(); >> 272 G4double x0den = material->GetIonisation()->GetX0density(); >> 273 G4double x1den = material->GetIonisation()->GetX1density(); 265 274 266 G4double eDensity = material->GetElectronDen 275 G4double eDensity = material->GetElectronDensity(); 267 276 268 // added ICRU90 stopping data for limited li << 277 G4double dedx = log(2.0*electron_mass_c2*bg2*cutEnergy/eexc2) 269 /* << 278 - (1.0 + cutEnergy/tmax)*beta2; 270 G4cout << "### DEDX ICRI90:" << (nullptr != << 271 << " Ekin=" << kineticEnergy << 272 << " " << p->GetParticleName() << 273 << " q2=" << chargeSquare << 274 << " inside " << material->GetName() << G4 << 275 */ << 276 if(nullptr != fICRU90 && kineticEnergy < fPr << 277 if(material != currentMaterial) { << 278 currentMaterial = material; << 279 baseMaterial = material->GetBaseMaterial << 280 ? material->GetBaseMaterial() : materi << 281 iICRU90 = fICRU90->GetIndex(baseMaterial << 282 } << 283 if(iICRU90 >= 0) { << 284 G4double dedx = 0.0; << 285 // only for alpha << 286 if(isAlpha) { << 287 if(kineticEnergy <= fAlphaTlimit) { << 288 dedx = fICRU90->GetElectronicDEDXforAlpha( << 289 } else { << 290 const G4double e = kineticEnergy*CLH << 291 dedx = fICRU90->GetElectronicDEDXforProton << 292 } << 293 } else { << 294 dedx = fICRU90->GetElectronicDEDXforPr << 295 *chargeSquare; << 296 } << 297 dedx *= material->GetDensity(); << 298 if(cutEnergy < tmax) { << 299 dedx += (G4Log(xc) + (1.0 - xc)*beta2) << 300 *(eDensity*chargeSquare/beta2); << 301 } << 302 //G4cout << " iICRU90=" << iICRU90 << << 303 if(dedx > 0.0) { return dedx; } << 304 } << 305 } << 306 // general Bethe-Bloch formula << 307 G4double dedx = G4Log(2.0*CLHEP::electron_ma << 308 - (1.0 + xc)*beta2; << 309 279 310 if(0.0 < spin) { << 280 if(0.5 == spin) { 311 G4double del = 0.5*cutEnergy/(kineticEnerg 281 G4double del = 0.5*cutEnergy/(kineticEnergy + mass); 312 dedx += del*del; 282 dedx += del*del; 313 } 283 } 314 284 315 // density correction 285 // density correction 316 G4double x = G4Log(bg2)/twoln10; << 286 G4double x = log(bg2)/twoln10; 317 dedx -= material->GetIonisation()->DensityCo << 287 if ( x >= x0den ) { >> 288 dedx -= twoln10*x - cden ; >> 289 if ( x < x1den ) dedx -= aden*pow((x1den-x),mden) ; >> 290 } 318 291 319 // shell correction 292 // shell correction 320 dedx -= 2.0*corr->ShellCorrection(p,material 293 dedx -= 2.0*corr->ShellCorrection(p,material,kineticEnergy); 321 294 322 // now compute the total ionization loss 295 // now compute the total ionization loss 323 dedx *= CLHEP::twopi_mc2_rcl2*chargeSquare*e << 296 >> 297 if (dedx < 0.0) dedx = 0.0 ; >> 298 >> 299 dedx *= twopi_mc2_rcl2*chargeSquare*eDensity/beta2; 324 300 325 //High order correction different for hadron 301 //High order correction different for hadrons and ions 326 if(isIon) { 302 if(isIon) { 327 dedx += corr->IonBarkasCorrection(p,materi 303 dedx += corr->IonBarkasCorrection(p,material,kineticEnergy); 328 } else { 304 } else { 329 dedx += corr->HighOrderCorrections(p,mater 305 dedx += corr->HighOrderCorrections(p,material,kineticEnergy,cutEnergy); 330 } 306 } 331 << 332 dedx = std::max(dedx, 0.0); << 333 /* << 334 G4cout << "E(MeV)= " << kineticEnergy/CLHEP: << 335 << " " << material->GetName() << G << 336 */ << 337 return dedx; 307 return dedx; 338 } 308 } 339 309 340 //....oooOO0OOooo........oooOO0OOooo........oo 310 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 341 311 342 void G4BetheBlochModel::CorrectionsAlongStep(c 312 void G4BetheBlochModel::CorrectionsAlongStep(const G4MaterialCutsCouple* couple, 343 c << 313 const G4DynamicParticle* dp, 344 c << 314 G4double& eloss, 345 G << 315 G4double&, 346 { << 316 G4double length) 347 // no correction for alpha << 317 { 348 if(isAlpha) { return; } << 349 << 350 // no correction at the last step or at smal << 351 const G4double preKinEnergy = dp->GetKinetic << 352 if(eloss >= preKinEnergy || eloss < preKinEn << 353 << 354 // corrections for all charged particles wit << 355 const G4ParticleDefinition* p = dp->GetDefin 318 const G4ParticleDefinition* p = dp->GetDefinition(); 356 if(p != particle) { SetupParameters(p); } << 357 if(!isIon) { return; } << 358 << 359 // effective energy and charge at a step << 360 const G4double e = std::max(preKinEnergy - e << 361 const G4Material* mat = couple->GetMaterial( 319 const G4Material* mat = couple->GetMaterial(); 362 const G4double q20 = corr->EffectiveChargeSq << 320 G4double preKinEnergy = dp->GetKineticEnergy(); 363 const G4double q2 = corr->EffectiveChargeSqu << 321 G4double e = preKinEnergy - eloss*0.5; 364 const G4double qfactor = q2/q20; << 322 if(e < 0.0) e = preKinEnergy*0.5; 365 << 323 366 /* << 324 if(isIon) { 367 G4cout << "G4BetheBlochModel::CorrectionsA << 325 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,e); 368 << preKinEnergy << " Eeff(MeV)=" << e << 326 GetModelOfFluctuations()->SetParticleAndCharge(p, q2); 369 << " eloss=" << eloss << " elossnew=" << e << 327 eloss *= q2*corr->EffectiveChargeCorrection(p,mat,e)/corrFactor; 370 << " qfactor=" << qfactor << " Qpre=" << q << 328 eloss += length*corr->IonHighOrderCorrections(p,couple,e); 371 << p->GetParticleName() <<G4endl; << 329 } 372 */ << 330 373 eloss *= qfactor; << 331 if(nuclearStopping && preKinEnergy*proton_mass_c2/mass < chargeSquare*100.*MeV) { >> 332 >> 333 G4double nloss = length*corr->NuclearDEDX(p,mat,e,false); >> 334 >> 335 // too big energy loss >> 336 if(eloss + nloss > preKinEnergy) { >> 337 nloss *= (preKinEnergy/(eloss + nloss)); >> 338 eloss = preKinEnergy; >> 339 } else { >> 340 eloss += nloss; >> 341 } >> 342 /* >> 343 G4cout << "G4ionIonisation::CorrectionsAlongStep: e= " << preKinEnergy >> 344 << " de= " << eloss << " NIEL= " << nloss >> 345 << " dynQ= " << dp->GetCharge()/eplus << G4endl; >> 346 */ >> 347 fParticleChange->ProposeNonIonizingEnergyDeposit(nloss); >> 348 } >> 349 374 } 350 } 375 351 376 //....oooOO0OOooo........oooOO0OOooo........oo 352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 377 353 378 void G4BetheBlochModel::SampleSecondaries(std: << 354 void G4BetheBlochModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp, 379 cons << 355 const G4MaterialCutsCouple*, 380 cons << 356 const G4DynamicParticle* dp, 381 G4do << 357 G4double minKinEnergy, 382 G4do << 358 G4double maxEnergy) 383 { << 359 { 384 G4double kinEnergy = dp->GetKineticEnergy(); << 360 G4double kineticEnergy = dp->GetKineticEnergy(); 385 const G4double tmax = MaxSecondaryEnergy(dp- << 361 G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(),kineticEnergy); 386 const G4double minKinEnergy = std::min(cut, << 362 387 const G4double maxKinEnergy = std::min(maxEn << 363 G4double maxKinEnergy = std::min(maxEnergy,tmax); 388 if(minKinEnergy >= maxKinEnergy) { return; } << 364 if(minKinEnergy >= maxKinEnergy) return; 389 << 365 390 //G4cout << "G4BetheBlochModel::SampleSecond << 366 G4double totEnergy = kineticEnergy + mass; 391 // << " Emax= " << maxKinEnergy << G << 367 G4double etot2 = totEnergy*totEnergy; 392 << 368 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2; 393 const G4double totEnergy = kinEnergy + mass; << 394 const G4double etot2 = totEnergy*totEnergy; << 395 const G4double beta2 = kinEnergy*(kinEnergy << 396 369 397 G4double deltaKinEnergy, f; 370 G4double deltaKinEnergy, f; 398 G4double f1 = 0.0; 371 G4double f1 = 0.0; 399 G4double fmax = 1.0; 372 G4double fmax = 1.0; 400 if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy* << 373 if( 0.5 == spin ) fmax += 0.5*maxKinEnergy*maxKinEnergy/etot2; 401 << 402 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra << 403 G4double rndm[2]; << 404 374 405 // sampling without nuclear size effect 375 // sampling without nuclear size effect 406 do { 376 do { 407 rndmEngineMod->flatArray(2, rndm); << 377 G4double q = G4UniformRand(); 408 deltaKinEnergy = minKinEnergy*maxKinEnergy 378 deltaKinEnergy = minKinEnergy*maxKinEnergy 409 /(minKinEnergy*(1.0 - rndm << 379 /(minKinEnergy*(1.0 - q) + maxKinEnergy*q); 410 380 411 f = 1.0 - beta2*deltaKinEnergy/tmax; 381 f = 1.0 - beta2*deltaKinEnergy/tmax; 412 if( 0.0 < spin ) { << 382 if( 0.5 == spin ) { 413 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/e 383 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2; 414 f += f1; 384 f += f1; 415 } 385 } 416 386 417 // Loop checking, 03-Aug-2015, Vladimir Iv << 387 } while( fmax*G4UniformRand() > f); 418 } while( fmax*rndm[1] > f); << 419 388 420 // projectile formfactor - suppresion of hig 389 // projectile formfactor - suppresion of high energy 421 // delta-electron production at high energy 390 // delta-electron production at high energy 422 391 423 G4double x = formfact*deltaKinEnergy; 392 G4double x = formfact*deltaKinEnergy; 424 if(x > 1.e-6) { 393 if(x > 1.e-6) { 425 394 426 G4double x1 = 1.0 + x; 395 G4double x1 = 1.0 + x; 427 G4double grej = 1.0/(x1*x1); << 396 G4double g = 1.0/(x1*x1); 428 if( 0.0 < spin ) { << 397 if( 0.5 == spin ) { 429 G4double x2 = 0.5*electron_mass_c2*delta 398 G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass); 430 grej *= (1.0 + magMoment2*(x2 - f1/f)/(1 << 399 g *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2)); 431 } 400 } 432 if(grej > 1.1) { << 401 if(g > 1.0) { 433 G4cout << "### G4BetheBlochModel WARNING << 402 G4cout << "### G4BetheBlochModel WARNING: g= " << g 434 << " " << dp->GetDefinition()->G << 403 << dp->GetDefinition()->GetParticleName() 435 << " Ekin(MeV)= " << kinEnergy << 404 << " Ekin(MeV)= " << kineticEnergy 436 << " delEkin(MeV)= " << deltaKinE << 405 << " delEkin(MeV)= " << deltaKinEnergy 437 << G4endl; << 406 << G4endl; 438 } 407 } 439 if(rndmEngineMod->flat() > grej) { return; << 408 if(G4UniformRand() > g) return; 440 } 409 } 441 410 442 G4ThreeVector deltaDirection; << 411 // delta-electron is produced 443 << 412 G4double totMomentum = totEnergy*sqrt(beta2); 444 if(UseAngularGeneratorFlag()) { << 413 G4double deltaMomentum = 445 const G4Material* mat = couple->GetMateria << 414 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2)); 446 deltaDirection = << 415 G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) / 447 GetAngularDistribution()->SampleDirectio << 416 (deltaMomentum * totMomentum); 448 SelectRandomAtomNumber(mat), << 417 if(cost > 1.0) { 449 mat); << 418 G4cout << "### G4BetheBlochModel WARNING: cost= " 450 } else { << 419 << cost << " > 1 for " 451 << 420 << dp->GetDefinition()->GetParticleName() 452 G4double deltaMomentum = << 421 << " Ekin(MeV)= " << kineticEnergy 453 std::sqrt(deltaKinEnergy * (deltaKinEner << 422 << " p(MeV/c)= " << totMomentum 454 G4double cost = deltaKinEnergy * (totEnerg << 423 << " delEkin(MeV)= " << deltaKinEnergy 455 (deltaMomentum * dp->GetTotalMomentum()) << 424 << " delMom(MeV/c)= " << deltaMomentum 456 cost = std::min(cost, 1.0); << 425 << " tmin(MeV)= " << minKinEnergy 457 const G4double sint = std::sqrt((1.0 - cos << 426 << " tmax(MeV)= " << maxKinEnergy 458 const G4double phi = twopi*rndmEngineMod-> << 459 << 460 deltaDirection.set(sint*std::cos(phi),sint << 461 deltaDirection.rotateUz(dp->GetMomentumDir << 462 } << 463 /* << 464 G4cout << "### G4BetheBlochModel " << 465 << dp->GetDefinition()->GetParticle << 466 << " Ekin(MeV)= " << kinEnergy << 467 << " delEkin(MeV)= " << deltaKinEne << 468 << " tmin(MeV)= " << minKinEnergy << 469 << " tmax(MeV)= " << maxKinEnergy << 470 << " dir= " << dp->GetMomentumDirec 427 << " dir= " << dp->GetMomentumDirection() 471 << " dirDelta= " << deltaDirection << 428 << G4endl; 472 << G4endl; << 429 cost = 1.0; 473 */ << 430 } >> 431 G4double sint = sqrt((1.0 - cost)*(1.0 + cost)); >> 432 >> 433 G4double phi = twopi * G4UniformRand() ; >> 434 >> 435 >> 436 G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost); >> 437 G4ThreeVector direction = dp->GetMomentumDirection(); >> 438 deltaDirection.rotateUz(direction); >> 439 474 // create G4DynamicParticle object for delta 440 // create G4DynamicParticle object for delta ray 475 auto delta = new G4DynamicParticle(theElectr << 441 G4DynamicParticle* delta = new G4DynamicParticle(theElectron, >> 442 deltaDirection,deltaKinEnergy); 476 443 477 vdp->push_back(delta); 444 vdp->push_back(delta); 478 445 479 // Change kinematics of primary particle 446 // Change kinematics of primary particle 480 kinEnergy -= deltaKinEnergy; << 447 kineticEnergy -= deltaKinEnergy; 481 G4ThreeVector finalP = dp->GetMomentum() - d << 448 G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum; 482 finalP = finalP.unit(); << 449 finalP = finalP.unit(); 483 450 484 fParticleChange->SetProposedKineticEnergy(ki << 451 fParticleChange->SetProposedKineticEnergy(kineticEnergy); 485 fParticleChange->SetProposedMomentumDirectio 452 fParticleChange->SetProposedMomentumDirection(finalP); 486 } << 487 << 488 //....oooOO0OOooo........oooOO0OOooo........oo << 489 << 490 G4double G4BetheBlochModel::MaxSecondaryEnergy << 491 << 492 { << 493 // here particle type is checked for the cas << 494 // when this model is shared between particl << 495 if(pd != particle) { SetupParameters(pd); } << 496 G4double tau = kinEnergy/mass; << 497 return 2.0*CLHEP::electron_mass_c2*tau*(tau << 498 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio); << 499 } 453 } 500 454 501 //....oooOO0OOooo........oooOO0OOooo........oo 455 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 502 456