Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // ------------------------------------------- 27 // 28 // GEANT4 Class header file 29 // 30 // 31 // File name: G4BetheBlochModel 32 // 33 // Author: Vladimir Ivanchenko on base 34 // 35 // Creation date: 03.01.2002 36 // 37 // Modifications: 38 // 39 // 04-12-02 Fix problem of G4DynamicParticle c 40 // 23-12-02 Change interface in order to move 41 // 27-01-03 Make models region aware (V.Ivanch 42 // 13-02-03 Add name (V.Ivanchenko) 43 // 24-03-05 Add G4EmCorrections (V.Ivanchenko) 44 // 11-04-05 Major optimisation of internal int 45 // 11-02-06 ComputeCrossSectionPerElectron, Co 46 // 12-02-06 move G4LossTableManager::Instance( 47 // in constructor (mma) 48 // 12-08-08 Added methods GetParticleCharge, G 49 // CorrectionsAlongStep needed for io 50 // 51 // ------------------------------------------- 52 // 53 54 //....oooOO0OOooo........oooOO0OOooo........oo 55 //....oooOO0OOooo........oooOO0OOooo........oo 56 57 #include "G4BetheBlochModel.hh" 58 #include "Randomize.hh" 59 #include "G4PhysicalConstants.hh" 60 #include "G4SystemOfUnits.hh" 61 #include "G4NistManager.hh" 62 #include "G4Electron.hh" 63 #include "G4LossTableManager.hh" 64 #include "G4EmCorrections.hh" 65 #include "G4EmParameters.hh" 66 #include "G4ParticleChangeForLoss.hh" 67 #include "G4ICRU90StoppingData.hh" 68 #include "G4Log.hh" 69 #include "G4DeltaAngle.hh" 70 #include <vector> 71 72 //....oooOO0OOooo........oooOO0OOooo........oo 73 74 G4BetheBlochModel::G4BetheBlochModel(const G4P 75 const G4S 76 : G4VEmModel(nam), 77 twoln10(2.0*G4Log(10.0)), 78 fAlphaTlimit(1*CLHEP::GeV), 79 fProtonTlimit(10*CLHEP::GeV) 80 { 81 theElectron = G4Electron::Electron(); 82 corr = G4LossTableManager::Instance()->EmCor 83 nist = G4NistManager::Instance(); 84 SetLowEnergyLimit(2.0*CLHEP::MeV); 85 } 86 87 //....oooOO0OOooo........oooOO0OOooo........oo 88 89 G4BetheBlochModel::~G4BetheBlochModel() = defa 90 91 //....oooOO0OOooo........oooOO0OOooo........oo 92 93 void G4BetheBlochModel::Initialise(const G4Par 94 const G4Dat 95 { 96 if(p != particle) { SetupParameters(p); } 97 98 // always false before the run 99 SetDeexcitationFlag(false); 100 101 // initialisation once 102 if(nullptr == fParticleChange) { 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 fParticleChange = GetParticleChangeForLoss 117 if(UseAngularGeneratorFlag() && nullptr == 118 SetAngularDistribution(new G4DeltaAngle( 119 } 120 } 121 // initialisation for each new run 122 if(IsMaster() && nullptr != fICRU90) { 123 fICRU90->Initialise(); 124 } 125 } 126 127 //....oooOO0OOooo........oooOO0OOooo........oo 128 129 G4double G4BetheBlochModel::GetChargeSquareRat 130 131 132 { 133 // this method is called only for ions, so n 134 if(isAlpha) { return 1.0; } 135 chargeSquare = corr->EffectiveChargeSquareRa 136 return chargeSquare; 137 } 138 139 //....oooOO0OOooo........oooOO0OOooo........oo 140 141 G4double G4BetheBlochModel::GetParticleCharge( 142 143 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 { 153 particle = p; 154 mass = particle->GetPDGMass(); 155 spin = particle->GetPDGSpin(); 156 G4double q = particle->GetPDGCharge()*invepl 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 } 175 176 //....oooOO0OOooo........oooOO0OOooo........oo 177 178 G4double G4BetheBlochModel::MinEnergyCut(const 179 const 180 { 181 return couple->GetMaterial()->GetIonisation( 182 } 183 184 //....oooOO0OOooo........oooOO0OOooo........oo 185 186 G4double 187 G4BetheBlochModel::ComputeCrossSectionPerElect 188 189 190 191 { 192 G4double cross = 0.0; 193 const G4double tmax = MaxSecondaryEnergy(p, 194 const G4double cutEnergy = std::min(std::min 195 const G4double maxEnergy = std::min(tmax, ma 196 if(cutEnergy < maxEnergy) { 197 198 G4double totEnergy = kineticEnergy + mass; 199 G4double energy2 = totEnergy*totEnergy; 200 G4double beta2 = kineticEnergy*(kineti 201 202 cross = (maxEnergy - cutEnergy)/(cutEnergy 203 - beta2*G4Log(maxEnergy/cutEnergy)/tmax; 204 205 // +term for spin=1/2 particle 206 if( 0.0 < spin ) { cross += 0.5*(maxEnergy 207 208 cross *= CLHEP::twopi_mc2_rcl2*chargeSquar 209 } 210 211 // G4cout << "BB: e= " << kineticEnergy << 212 // << " tmax= " << tmax << " cross= 213 214 return cross; 215 } 216 217 //....oooOO0OOooo........oooOO0OOooo........oo 218 219 G4double G4BetheBlochModel::ComputeCrossSectio 220 con 221 222 223 224 225 { 226 return Z*ComputeCrossSectionPerElectron(p,ki 227 } 228 229 //....oooOO0OOooo........oooOO0OOooo........oo 230 231 G4double G4BetheBlochModel::CrossSectionPerVol 232 con 233 con 234 235 236 237 { 238 G4double sigma = mat->GetElectronDensity() 239 *ComputeCrossSectionPerElectron(p,kinEnerg 240 if(isAlpha) { 241 sigma *= corr->EffectiveChargeSquareRatio( 242 } 243 return sigma; 244 } 245 246 //....oooOO0OOooo........oooOO0OOooo........oo 247 248 G4double G4BetheBlochModel::ComputeDEDXPerVolu 249 250 251 252 { 253 const G4double tmax = MaxSecondaryEnergy(p, 254 // projectile formfactor limit energy loss 255 const G4double cutEnergy = std::min(std::min 256 257 G4double tau = kineticEnergy/mass; 258 G4double gam = tau + 1.0; 259 G4double bg2 = tau * (tau+2.0); 260 G4double beta2 = bg2/(gam*gam); 261 G4double xc = cutEnergy/tmax; 262 263 G4double eexc = material->GetIonisation()-> 264 G4double eexc2 = eexc*eexc; 265 266 G4double eDensity = material->GetElectronDen 267 268 // added ICRU90 stopping data for limited li 269 /* 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 310 if(0.0 < spin) { 311 G4double del = 0.5*cutEnergy/(kineticEnerg 312 dedx += del*del; 313 } 314 315 // density correction 316 G4double x = G4Log(bg2)/twoln10; 317 dedx -= material->GetIonisation()->DensityCo 318 319 // shell correction 320 dedx -= 2.0*corr->ShellCorrection(p,material 321 322 // now compute the total ionization loss 323 dedx *= CLHEP::twopi_mc2_rcl2*chargeSquare*e 324 325 //High order correction different for hadron 326 if(isIon) { 327 dedx += corr->IonBarkasCorrection(p,materi 328 } else { 329 dedx += corr->HighOrderCorrections(p,mater 330 } 331 332 dedx = std::max(dedx, 0.0); 333 /* 334 G4cout << "E(MeV)= " << kineticEnergy/CLHEP: 335 << " " << material->GetName() << G 336 */ 337 return dedx; 338 } 339 340 //....oooOO0OOooo........oooOO0OOooo........oo 341 342 void G4BetheBlochModel::CorrectionsAlongStep(c 343 c 344 c 345 G 346 { 347 // no correction for alpha 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 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( 362 const G4double q20 = corr->EffectiveChargeSq 363 const G4double q2 = corr->EffectiveChargeSqu 364 const G4double qfactor = q2/q20; 365 366 /* 367 G4cout << "G4BetheBlochModel::CorrectionsA 368 << preKinEnergy << " Eeff(MeV)=" << e 369 << " eloss=" << eloss << " elossnew=" << e 370 << " qfactor=" << qfactor << " Qpre=" << q 371 << p->GetParticleName() <<G4endl; 372 */ 373 eloss *= qfactor; 374 } 375 376 //....oooOO0OOooo........oooOO0OOooo........oo 377 378 void G4BetheBlochModel::SampleSecondaries(std: 379 cons 380 cons 381 G4do 382 G4do 383 { 384 G4double kinEnergy = dp->GetKineticEnergy(); 385 const G4double tmax = MaxSecondaryEnergy(dp- 386 const G4double minKinEnergy = std::min(cut, 387 const G4double maxKinEnergy = std::min(maxEn 388 if(minKinEnergy >= maxKinEnergy) { return; } 389 390 //G4cout << "G4BetheBlochModel::SampleSecond 391 // << " Emax= " << maxKinEnergy << G 392 393 const G4double totEnergy = kinEnergy + mass; 394 const G4double etot2 = totEnergy*totEnergy; 395 const G4double beta2 = kinEnergy*(kinEnergy 396 397 G4double deltaKinEnergy, f; 398 G4double f1 = 0.0; 399 G4double fmax = 1.0; 400 if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy* 401 402 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra 403 G4double rndm[2]; 404 405 // sampling without nuclear size effect 406 do { 407 rndmEngineMod->flatArray(2, rndm); 408 deltaKinEnergy = minKinEnergy*maxKinEnergy 409 /(minKinEnergy*(1.0 - rndm 410 411 f = 1.0 - beta2*deltaKinEnergy/tmax; 412 if( 0.0 < spin ) { 413 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/e 414 f += f1; 415 } 416 417 // Loop checking, 03-Aug-2015, Vladimir Iv 418 } while( fmax*rndm[1] > f); 419 420 // projectile formfactor - suppresion of hig 421 // delta-electron production at high energy 422 423 G4double x = formfact*deltaKinEnergy; 424 if(x > 1.e-6) { 425 426 G4double x1 = 1.0 + x; 427 G4double grej = 1.0/(x1*x1); 428 if( 0.0 < spin ) { 429 G4double x2 = 0.5*electron_mass_c2*delta 430 grej *= (1.0 + magMoment2*(x2 - f1/f)/(1 431 } 432 if(grej > 1.1) { 433 G4cout << "### G4BetheBlochModel WARNING 434 << " " << dp->GetDefinition()->G 435 << " Ekin(MeV)= " << kinEnergy 436 << " delEkin(MeV)= " << deltaKinE 437 << G4endl; 438 } 439 if(rndmEngineMod->flat() > grej) { return; 440 } 441 442 G4ThreeVector deltaDirection; 443 444 if(UseAngularGeneratorFlag()) { 445 const G4Material* mat = couple->GetMateria 446 deltaDirection = 447 GetAngularDistribution()->SampleDirectio 448 SelectRandomAtomNumber(mat), 449 mat); 450 } else { 451 452 G4double deltaMomentum = 453 std::sqrt(deltaKinEnergy * (deltaKinEner 454 G4double cost = deltaKinEnergy * (totEnerg 455 (deltaMomentum * dp->GetTotalMomentum()) 456 cost = std::min(cost, 1.0); 457 const G4double sint = std::sqrt((1.0 - cos 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 471 << " dirDelta= " << deltaDirection 472 << G4endl; 473 */ 474 // create G4DynamicParticle object for delta 475 auto delta = new G4DynamicParticle(theElectr 476 477 vdp->push_back(delta); 478 479 // Change kinematics of primary particle 480 kinEnergy -= deltaKinEnergy; 481 G4ThreeVector finalP = dp->GetMomentum() - d 482 finalP = finalP.unit(); 483 484 fParticleChange->SetProposedKineticEnergy(ki 485 fParticleChange->SetProposedMomentumDirectio 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 } 500 501 //....oooOO0OOooo........oooOO0OOooo........oo 502