Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // ------------------------------------------------------------------- 27 // 28 // GEANT4 Class file 29 // 30 // 31 // File name: G4SeltzerBergerModel 32 // 33 // Author: Vladimir Ivanchenko use inheritance from Andreas Schaelicke 34 // base class implementing ultra relativistic bremsstrahlung 35 // model 36 // 37 // Creation date: 04.10.2011 38 // 39 // Modifications: 40 // 41 // 24.07.2018 Introduced possibility to use sampling tables to sample the 42 // emitted photon energy (instead of using rejectio) from the 43 // Seltzer-Berger scalled DCS for bremsstrahlung photon emission. 44 // Using these sampling tables option gives faster(30-70%) final 45 // state generation than the original rejection but takes some 46 // extra memory (+ ~6MB in the case of the full CMS detector). 47 // (M Novak) 48 // 49 // ------------------------------------------------------------------- 50 // 51 52 #include "G4SeltzerBergerModel.hh" 53 #include "G4PhysicalConstants.hh" 54 #include "G4SystemOfUnits.hh" 55 #include "Randomize.hh" 56 #include "G4Material.hh" 57 #include "G4Element.hh" 58 #include "G4ElementVector.hh" 59 #include "G4ParticleChangeForLoss.hh" 60 #include "G4SBBremTable.hh" 61 #include "G4ModifiedTsai.hh" 62 63 #include "G4EmParameters.hh" 64 #include "G4ProductionCutsTable.hh" 65 #include "G4Gamma.hh" 66 #include "G4Electron.hh" 67 68 #include "G4Physics2DVector.hh" 69 #include "G4Exp.hh" 70 #include "G4Log.hh" 71 #include "G4AutoLock.hh" 72 73 #include "G4ios.hh" 74 75 #include <fstream> 76 #include <iomanip> 77 #include <sstream> 78 #include <thread> 79 80 G4double G4SeltzerBergerModel::gYLimitData[] = { 0.0 }; 81 G4Physics2DVector* G4SeltzerBergerModel::gSBDCSData[] = { nullptr }; 82 G4SBBremTable* G4SeltzerBergerModel::gSBSamplingTable = nullptr; 83 84 // constant DCS factor: 16\alpha r_0^2/3 85 const G4double G4SeltzerBergerModel::gBremFactor 86 = 16. * CLHEP::fine_structure_const * CLHEP::classic_electr_radius 87 * CLHEP::classic_electr_radius/3.; 88 89 // Migdal's constant: 4\pi r_0*electron_reduced_compton_wavelength^2 90 const G4double G4SeltzerBergerModel::gMigdalConstant 91 = 4. * CLHEP::pi * CLHEP::classic_electr_radius 92 * CLHEP::electron_Compton_length * CLHEP::electron_Compton_length; 93 94 static constexpr G4double twoMass = 2* CLHEP::electron_mass_c2; 95 static constexpr G4double kAlpha = CLHEP::twopi*CLHEP::fine_structure_const; 96 static std::once_flag applyOnce; 97 98 namespace 99 { 100 G4Mutex theSBMutex = G4MUTEX_INITIALIZER; 101 102 // for numerical integration on [0,1] 103 const G4double gXGL[8] = { 104 1.98550718e-02, 1.01666761e-01, 2.37233795e-01, 4.08282679e-01, 105 5.91717321e-01, 7.62766205e-01, 8.98333239e-01, 9.80144928e-01 106 }; 107 const G4double gWGL[8] = { 108 5.06142681e-02, 1.11190517e-01, 1.56853323e-01, 1.81341892e-01, 109 1.81341892e-01, 1.56853323e-01, 1.11190517e-01, 5.06142681e-02 110 }; 111 } 112 113 G4SeltzerBergerModel::G4SeltzerBergerModel(const G4ParticleDefinition* p, 114 const G4String& nam) 115 : G4VEmModel(nam), 116 fGammaParticle(G4Gamma::Gamma()), 117 fLowestKinEnergy(1.0*CLHEP::keV) 118 { 119 SetLowEnergyLimit(fLowestKinEnergy); 120 SetAngularDistribution(new G4ModifiedTsai()); 121 if (fPrimaryParticle != p) { SetParticle(p); } 122 } 123 124 G4SeltzerBergerModel::~G4SeltzerBergerModel() 125 { 126 // delete SB-DCS data per Z 127 if (isInitializer) { 128 for (std::size_t iz = 0; iz < gMaxZet; ++iz) { 129 if (gSBDCSData[iz]) { 130 delete gSBDCSData[iz]; 131 gSBDCSData[iz] = nullptr; 132 } 133 } 134 if (gSBSamplingTable) { 135 delete gSBSamplingTable; 136 gSBSamplingTable = nullptr; 137 } 138 } 139 } 140 141 void G4SeltzerBergerModel::Initialise(const G4ParticleDefinition* p, 142 const G4DataVector& cuts) 143 { 144 // parameters in each thread 145 if (fPrimaryParticle != p) { 146 SetParticle(p); 147 } 148 fIsUseSamplingTables = G4EmParameters::Instance()->EnableSamplingTable(); 149 fCurrentIZ = 0; 150 151 // initialise static tables for the Seltzer-Berger model 152 std::call_once(applyOnce, [this]() { isInitializer = true; }); 153 154 if (isInitializer) { 155 G4AutoLock l(&theSBMutex); 156 157 // initialisation per element is done only once 158 auto elemTable = G4Element::GetElementTable(); 159 for (auto const & elm : *elemTable) { 160 G4int Z = std::max(1,std::min(elm->GetZasInt(), gMaxZet-1)); 161 // load SB-DCS data for this atomic number if it has not been loaded yet 162 if (gSBDCSData[Z] == nullptr) ReadData(Z); 163 } 164 165 // init sampling tables if it was requested 166 if (fIsUseSamplingTables) { 167 if (nullptr == gSBSamplingTable) { 168 gSBSamplingTable = new G4SBBremTable(); 169 } 170 gSBSamplingTable->Initialize(std::max(fLowestKinEnergy, LowEnergyLimit()), 171 HighEnergyLimit()); 172 } 173 l.unlock(); 174 } 175 // element selectors are initialized in the master thread 176 if (IsMaster()) { 177 InitialiseElementSelectors(p, cuts); 178 } 179 // initialisation in all threads 180 if (nullptr == fParticleChange) { 181 fParticleChange = GetParticleChangeForLoss(); 182 } 183 auto trmodel = GetTripletModel(); 184 if (nullptr != trmodel) { 185 trmodel->Initialise(p, cuts); 186 fIsScatOffElectron = true; 187 } 188 } 189 190 void G4SeltzerBergerModel::InitialiseLocal(const G4ParticleDefinition*, 191 G4VEmModel* masterModel) 192 { 193 SetElementSelectors(masterModel->GetElementSelectors()); 194 } 195 196 void G4SeltzerBergerModel::SetParticle(const G4ParticleDefinition* p) 197 { 198 fPrimaryParticle = p; 199 fIsElectron = (p == G4Electron::Electron()); 200 } 201 202 void G4SeltzerBergerModel::ReadData(G4int Z) { 203 // return if it has been already loaded 204 if (gSBDCSData[Z] != nullptr) return; 205 206 if (gSBDCSData[Z] == nullptr) { 207 std::ostringstream ost; 208 ost << G4EmParameters::Instance()->GetDirLEDATA() << "/brem_SB/br" << Z; 209 std::ifstream fin(ost.str().c_str()); 210 if (!fin.is_open()) { 211 G4ExceptionDescription ed; 212 ed << "Bremsstrahlung data file <" << ost.str().c_str() 213 << "> is not opened!"; 214 G4Exception("G4SeltzerBergerModel::ReadData()","em0003",FatalException, 215 ed,"G4LEDATA version should be G4EMLOW6.23 or later."); 216 return; 217 } 218 //G4cout << "G4SeltzerBergerModel read from <" << ost.str().c_str() 219 // << ">" << G4endl; 220 auto v = new G4Physics2DVector(); 221 if (v->Retrieve(fin)) { 222 v->SetBicubicInterpolation(fIsUseBicubicInterpolation); 223 static const G4double emaxlog = 4*G4Log(10.); 224 gYLimitData[Z] = v->Value(0.97, emaxlog, fIndx, fIndy); 225 gSBDCSData[Z] = v; 226 } else { 227 G4ExceptionDescription ed; 228 ed << "Bremsstrahlung data file <" << ost.str().c_str() 229 << "> is not retrieved!"; 230 G4Exception("G4SeltzerBergerModel::ReadData()","em0005",FatalException, 231 ed,"G4LEDATA version should be G4EMLOW6.23 or later."); 232 delete v; 233 } 234 } 235 } 236 237 // minimum primary (e-/e+) energy at which discrete interaction is possible 238 G4double G4SeltzerBergerModel::MinPrimaryEnergy(const G4Material*, 239 const G4ParticleDefinition*, 240 G4double cut) 241 { 242 return std::max(fLowestKinEnergy, cut); 243 } 244 245 // Sets kinematical variables like E_kin, E_t and some material dependent 246 // for characteristic photon energy k_p (more exactly 247 // k_p^2) for the Ter-Mikaelian suppression effect. 248 void G4SeltzerBergerModel::SetupForMaterial(const G4ParticleDefinition*, 249 const G4Material* mat, 250 G4double kinEnergy) 251 { 252 fDensityFactor = gMigdalConstant*mat->GetElectronDensity(); 253 // calculate threshold for density effect: k_p = sqrt(fDensityCorr) 254 fPrimaryKinEnergy = kinEnergy; 255 fPrimaryTotalEnergy = kinEnergy + CLHEP::electron_mass_c2; 256 fDensityCorr = fDensityFactor*fPrimaryTotalEnergy*fPrimaryTotalEnergy; 257 } 258 259 // Computes the restricted dE/dx as the appropriate weight of the individual 260 // element contributions that are computed by numerically integrating the DCS. 261 G4double 262 G4SeltzerBergerModel::ComputeDEDXPerVolume(const G4Material* material, 263 const G4ParticleDefinition* p, 264 G4double kineticEnergy, 265 G4double cutEnergy) 266 { 267 G4double dedx = 0.0; 268 if (nullptr == fPrimaryParticle) { 269 SetParticle(p); 270 } 271 if (kineticEnergy <= fLowestKinEnergy) { 272 return dedx; 273 } 274 // maximum value of the dE/dx integral (the minimum is 0 of course) 275 G4double tmax = std::min(cutEnergy, kineticEnergy); 276 if (tmax == 0.0) { 277 return dedx; 278 } 279 // sets kinematical and material related variables 280 SetupForMaterial(fPrimaryParticle, material, kineticEnergy); 281 // get element compositions of the material 282 const G4ElementVector* theElemVector = material->GetElementVector(); 283 const G4double* theAtomNumDensVector = material->GetAtomicNumDensityVector(); 284 const std::size_t numberOfElements = theElemVector->size(); 285 // loop over the elements of the material and compute their contributions to 286 // the restricted dE/dx by numerical integration of the dependent part of DCS 287 for (std::size_t ie = 0; ie < numberOfElements; ++ie) { 288 G4VEmModel::SetCurrentElement((*theElemVector)[ie]); 289 G4int Z = (*theElemVector)[ie]->GetZasInt(); 290 fCurrentIZ = std::min(Z, gMaxZet); 291 dedx += (Z*Z)*theAtomNumDensVector[ie]*ComputeBremLoss(tmax); 292 } 293 // apply the constant factor C/Z = 16\alpha r_0^2/3 294 dedx *= gBremFactor; 295 return std::max(dedx, 0.); 296 } 297 298 // Computes the integral part of the restricted dE/dx contribution from a given 299 // element (Z) by numerically integrating the k dependent DCS between 300 // k_min=0 and k_max = tmax = min[gamma-cut, electron-kinetic-energy]. 301 // The numerical integration is done by dividing the integration range into 'n' 302 // subintervals and an 8 pint GL integral (on [0,1]) is performed on each sub- 303 // inteval by tranforming k to alpha=k/E_t (E_t is the total energy of the e-) 304 // and each sub-interavl is transformed to [0,1]. So the integrastion is done 305 // in xi(alpha) = xi(k) = [k/E_t-alpha_i]/delta where alpha_i=(i-1)*delta for 306 // the i = 1,2,..,n-th sub-interval so xi(k) in [0,1] on each sub-intevals. 307 // This transformation from 'k' to 'xi(k)' results in a multiplicative factor 308 // of E_t*delta at each step. 309 // The restricted dE/dx = N int_{0}^{k_max} k*ds/dk dk. In this case not 310 // the ds/dk(Z,k) but ds/dk(Z,k)*[F*k/C] is computed since: 311 // (i) what we need here is ds/dk*k and not k so this multiplication is done 312 // (ii) the Ter-Mikaelian suppression i.e. F related factor is done here 313 // (iii) the constant factor C (includes Z^2 as well)is accounted in the caller 314 G4double G4SeltzerBergerModel::ComputeBremLoss(G4double tmax) 315 { 316 // number of intervals and integration step 317 const G4double alphaMax = tmax/fPrimaryTotalEnergy; 318 const G4int nSub = (G4int)(20*alphaMax)+3; 319 const G4double delta = alphaMax/((G4double)nSub); 320 // set minimum value of the first sub-inteval 321 G4double alpha_i = 0.0; 322 G4double dedxInteg = 0.0; 323 for (G4int l = 0; l < nSub; ++l) { 324 for (G4int igl = 0; igl < 8; ++igl) { 325 // compute the emitted photon energy k 326 const G4double k = (alpha_i+gXGL[igl]*delta)*fPrimaryTotalEnergy; 327 // compute the DCS value at k (without the constant, the 1/k, 1/F factors) 328 const G4double dcs = ComputeDXSectionPerAtom(k); 329 // account Ter-Mikaelian suppression: times 1/F with F = 1+(k_p/k)^2 330 dedxInteg += gWGL[igl]*dcs/(1.0+fDensityCorr/(k*k)); 331 } 332 // update sub-interval minimum value 333 alpha_i += delta; 334 } 335 // apply corrections due to variable transformation i.e. E_t*delta 336 dedxInteg *= delta*fPrimaryTotalEnergy; 337 return std::max(dedxInteg,0.); 338 } 339 340 // Computes restrected atomic cross section by numerically integrating the 341 // DCS between the proper kinematical limits accounting the gamma production cut 342 G4double 343 G4SeltzerBergerModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* p, 344 G4double kineticEnergy, 345 G4double Z, 346 G4double, 347 G4double cut, 348 G4double maxEnergy) 349 { 350 G4double crossSection = 0.0; 351 if (nullptr == fPrimaryParticle) { 352 SetParticle(p); 353 } 354 if (kineticEnergy <= fLowestKinEnergy) { 355 return crossSection; 356 } 357 // min/max kinetic energy limits of the DCS integration: 358 const G4double tmin = std::min(cut, kineticEnergy); 359 const G4double tmax = std::min(maxEnergy, kineticEnergy); 360 // zero restricted x-section if e- kinetic energy is below gamma cut 361 if (tmin >= tmax) { 362 return crossSection; 363 } 364 fCurrentIZ = std::min(G4lrint(Z), gMaxZet); 365 // integrate numerically (dependent part of) the DCS between the kin. limits: 366 // a. integrate between tmin and kineticEnergy of the e- 367 crossSection = ComputeXSectionPerAtom(tmin); 368 // allow partial integration: only if maxEnergy < kineticEnergy 369 // b. integrate between tmax and kineticEnergy (tmax=maxEnergy in this case) 370 // (so the result in this case is the integral of DCS between tmin and 371 // maxEnergy) 372 if (tmax < kineticEnergy) { 373 crossSection -= ComputeXSectionPerAtom(tmax); 374 } 375 // multiply with the constant factors: 16\alpha r_0^2/3 Z^2 376 crossSection *= Z*Z*gBremFactor; 377 return std::max(crossSection, 0.); 378 } 379 380 // Numerical integral of the (k dependent part of) DCS between k_min=tmin and 381 // k_max = E_k (where E_k is the kinetic energy of the e- and tmin is the 382 // minimum of energy of the emitted photon). The integration is done in the 383 // transformed alpha(k) = ln(k/E_t) variable (with E_t being the total energy of 384 // the primary e-). The integration range is divided into n sub-intervals with 385 // delta = [ln(k_min/E_t)-ln(k_max/E_t)]/n width each. An 8 point GL integral 386 // on [0,1] is applied on each sub-inteval so alpha is transformed to 387 // xi(alpha) = xi(k) = [ln(k/E_t)-alpha_i]/delta where alpha_i = ln(k_min/E_t) + 388 // (i-1)*delta for the i = 1,2,..,n-th sub-interval and xi(k) in [0,1] on each 389 // sub-intevals. From the transformed xi, k(xi) = E_t exp[xi*delta+alpha_i]. 390 // Since the integration is done in variable xi instead of k this 391 // transformation results in a multiplicative factor of k*delta at each step. 392 // However, DCS differential in k is ~1/k so the multiplicative factor is simple 393 // becomes delta and the 1/k factor is dropped from the DCS computation. 394 // Ter-Mikaelian suppression is always accounted 395 G4double G4SeltzerBergerModel::ComputeXSectionPerAtom(G4double tmin) 396 { 397 G4double xSection = 0.0; 398 const G4double alphaMin = G4Log(tmin/fPrimaryTotalEnergy); 399 const G4double alphaMax = G4Log(fPrimaryKinEnergy/fPrimaryTotalEnergy); 400 const G4int nSub = (G4int)(0.45*(alphaMax-alphaMin))+4; 401 const G4double delta = (alphaMax-alphaMin)/((G4double)nSub); 402 // set minimum value of the first sub-inteval 403 G4double alpha_i = alphaMin; 404 for (G4int l = 0; l < nSub; ++l) { 405 for (G4int igl = 0; igl < 8; ++igl) { 406 // compute the emitted photon energy k 407 const G4double k = G4Exp(alpha_i+gXGL[igl]*delta)*fPrimaryTotalEnergy; 408 // compute the DCS value at k (without the constant, the 1/k, 1/F factors) 409 const G4double dcs = ComputeDXSectionPerAtom(k); 410 // account Ter-Mikaelian suppression: times 1/F with F = 1+(k_p/k)^2 411 xSection += gWGL[igl]*dcs/(1.0+fDensityCorr/(k*k)); 412 } 413 // update sub-interval minimum value 414 alpha_i += delta; 415 } 416 // apply corrections due to variable transformation 417 xSection *= delta; 418 // final check 419 return std::max(xSection, 0.); 420 } 421 422 G4double G4SeltzerBergerModel::ComputeDXSectionPerAtom(G4double gammaEnergy) 423 { 424 G4double dxsec = 0.0; 425 if (gammaEnergy < 0.0 || fPrimaryKinEnergy <= 0.0) { 426 return dxsec; 427 } 428 // reduced photon energy 429 const G4double x = gammaEnergy/fPrimaryKinEnergy; 430 // l-kinetic energy of the e-/e+ 431 const G4double y = G4Log(fPrimaryKinEnergy/CLHEP::MeV); 432 // make sure that the Z-related SB-DCS are loaded 433 // NOTE: fCurrentIZ should have been set before. 434 fCurrentIZ = std::max(std::min(fCurrentIZ, gMaxZet-1), 1); 435 if (nullptr == gSBDCSData[fCurrentIZ]) { 436 G4AutoLock l(&theSBMutex); 437 ReadData(fCurrentIZ); 438 l.unlock(); 439 } 440 // NOTE: SetupForMaterial should have been called before! 441 const G4double pt2 = fPrimaryKinEnergy*(fPrimaryKinEnergy + twoMass); 442 const G4double invb2 = fPrimaryTotalEnergy*fPrimaryTotalEnergy/pt2; 443 G4double val = gSBDCSData[fCurrentIZ]->Value(x,y,fIndx,fIndy); 444 dxsec = val*invb2*CLHEP::millibarn/gBremFactor; 445 // e+ correction 446 if (!fIsElectron) { 447 const G4double invbeta1 = std::sqrt(invb2); 448 const G4double e2 = fPrimaryKinEnergy - gammaEnergy; 449 if (e2 > 0.0) { 450 const G4double invbeta2 = 451 (e2 + CLHEP::electron_mass_c2)/std::sqrt(e2*(e2 + twoMass)); 452 const G4double dum0 = kAlpha*fCurrentIZ*(invbeta1-invbeta2); 453 if (dum0 < gExpNumLimit) { 454 dxsec = 0.0; 455 } else { 456 dxsec *= G4Exp(dum0); 457 } 458 } else { 459 dxsec = 0.0; 460 } 461 } 462 return dxsec; 463 } 464 465 void 466 G4SeltzerBergerModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp, 467 const G4MaterialCutsCouple* couple, 468 const G4DynamicParticle* dp, 469 G4double cutEnergy, 470 G4double maxEnergy) 471 { 472 const G4double kinEnergy = dp->GetKineticEnergy(); 473 const G4double logKinEnergy = dp->GetLogKineticEnergy(); 474 const G4double tmin = std::min(cutEnergy, kinEnergy); 475 const G4double tmax = std::min(maxEnergy, kinEnergy); 476 if (tmin >= tmax) { 477 return; 478 } 479 // set local variables and select target element 480 SetupForMaterial(fPrimaryParticle, couple->GetMaterial(), kinEnergy); 481 const G4Element* elm = SelectTargetAtom(couple, fPrimaryParticle, kinEnergy, 482 logKinEnergy, tmin, tmax); 483 fCurrentIZ = std::max(std::min(elm->GetZasInt(), gMaxZet-1), 1); 484 // 485 const G4double totMomentum = std::sqrt(kinEnergy*(kinEnergy + twoMass)); 486 /* 487 G4cout << "G4SeltzerBergerModel::SampleSecondaries E(MeV)= " 488 << kinEnergy/MeV 489 << " Z= " << fCurrentIZ << " cut(MeV)= " << tmin/MeV 490 << " emax(MeV)= " << tmax/MeV << " corr= " << fDensityCorr << G4endl; 491 */ 492 // sample emitted photon energy either by rejection or from samplign tables 493 const G4double gammaEnergy = !fIsUseSamplingTables 494 ? SampleEnergyTransfer(kinEnergy, logKinEnergy, tmin, tmax) 495 : gSBSamplingTable->SampleEnergyTransfer(kinEnergy, logKinEnergy, tmin, 496 fDensityCorr, fCurrentIZ, couple->GetIndex(), fIsElectron); 497 // should never happen under normal conditions but protect it 498 if (gammaEnergy <= 0.) { 499 return; 500 } 501 // 502 // angles of the emitted gamma. ( Z - axis along the parent particle) use 503 // general interface 504 G4ThreeVector gamDir = GetAngularDistribution()->SampleDirection(dp, 505 fPrimaryTotalEnergy-gammaEnergy, fCurrentIZ, couple->GetMaterial()); 506 // create G4DynamicParticle object for the emitted Gamma 507 auto gamma = new G4DynamicParticle(fGammaParticle, gamDir, gammaEnergy); 508 vdp->push_back(gamma); 509 // 510 // compute post-interaction kinematics of the primary e-/e+ 511 G4ThreeVector dir = 512 (totMomentum*dp->GetMomentumDirection() - gammaEnergy*gamDir).unit(); 513 const G4double finalE = kinEnergy - gammaEnergy; 514 /* 515 G4cout << "### G4SBModel: v= " 516 << " Eg(MeV)= " << gammaEnergy 517 << " Ee(MeV)= " << kineticEnergy 518 << " DirE " << direction << " DirG " << gammaDirection 519 << G4endl; 520 */ 521 // if secondary gamma energy is higher than threshold(very high by default) 522 // then stop tracking the primary particle and create new secondary e-/e+ 523 // instead of the primary 524 if (gammaEnergy > SecondaryThreshold()) { 525 fParticleChange->ProposeTrackStatus(fStopAndKill); 526 fParticleChange->SetProposedKineticEnergy(0.0); 527 auto el = new G4DynamicParticle( 528 const_cast<G4ParticleDefinition*>(fPrimaryParticle), dir, finalE); 529 vdp->push_back(el); 530 } else { // continue tracking the primary e-/e+ otherwise 531 fParticleChange->SetProposedMomentumDirection(dir); 532 fParticleChange->SetProposedKineticEnergy(finalE); 533 } 534 } 535 536 // sample emitted photon energy by usign rejection 537 G4double 538 G4SeltzerBergerModel::SampleEnergyTransfer(const G4double kinEnergy, 539 const G4double logKinEnergy, 540 const G4double tmin, 541 const G4double tmax) 542 { 543 // min max of the transformed variable: x(k) = ln(k^2+k_p^2) that is in 544 // [ln(k_c^2+k_p^2), ln(E_k^2+k_p^2)] 545 const G4double xmin = G4Log(tmin*tmin+fDensityCorr); 546 const G4double xrange = G4Log(tmax*tmax+fDensityCorr)-xmin; 547 const G4double y = logKinEnergy; 548 // majoranta 549 const G4double x0 = tmin/kinEnergy; 550 G4double vmax; 551 if (nullptr == gSBDCSData[fCurrentIZ]) { 552 ReadData(fCurrentIZ); 553 } 554 vmax = gSBDCSData[fCurrentIZ]->Value(x0, y, fIndx, fIndy)*1.02; 555 // 556 static const G4double kEPeakLim = 300.*CLHEP::MeV; 557 static const G4double kELowLim = 20.*CLHEP::keV; 558 // majoranta corrected for e- 559 if (fIsElectron && x0 < 0.97 && 560 ((kinEnergy>kEPeakLim) || (kinEnergy<kELowLim))) { 561 G4double ylim = std::min(gYLimitData[fCurrentIZ], 562 1.1*gSBDCSData[fCurrentIZ]->Value(0.97,y,fIndx,fIndy)); 563 vmax = std::max(vmax, ylim); 564 } 565 if (x0 < 0.05) { 566 vmax *= 1.2; 567 } 568 //G4cout<<"y= "<<y<<" xmin= "<<xmin<<" xmax= "<<xmax 569 //<<" vmax= "<<vmax<<G4endl; 570 static const G4int kNCountMax = 100; 571 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine(); 572 G4double rndm[2]; 573 G4double gammaEnergy, v; 574 for (G4int nn = 0; nn < kNCountMax; ++nn) { 575 rndmEngine->flatArray(2, rndm); 576 gammaEnergy = 577 std::sqrt(std::max(G4Exp(xmin + rndm[0]*xrange)-fDensityCorr,0.)); 578 v = gSBDCSData[fCurrentIZ]->Value(gammaEnergy/kinEnergy, y, fIndx, fIndy); 579 // e+ correction 580 if (!fIsElectron) { 581 const G4double e1 = kinEnergy - tmin; 582 const G4double invbeta1 = 583 (e1 + CLHEP::electron_mass_c2)/std::sqrt(e1*(e1 + twoMass)); 584 const G4double e2 = kinEnergy-gammaEnergy; 585 const G4double invbeta2 = 586 (e2 + CLHEP::electron_mass_c2)/std::sqrt(e2*(e2 + twoMass)); 587 const G4double dum0 = kAlpha*fCurrentIZ*(invbeta1-invbeta2); 588 if (dum0 < gExpNumLimit) { 589 v = 0.0; 590 } else { 591 v *= G4Exp(dum0); 592 } 593 } 594 if (v > 1.05*vmax && fNumWarnings < 11) { 595 ++fNumWarnings; 596 G4ExceptionDescription ed; 597 ed << "### G4SeltzerBergerModel Warning: Majoranta exceeded! " 598 << v << " > " << vmax << " by " << v/vmax 599 << " Niter= " << nn 600 << " Egamma(MeV)= " << gammaEnergy 601 << " Ee(MeV)= " << kinEnergy 602 << " Z= " << fCurrentIZ << " " << fPrimaryParticle->GetParticleName(); 603 // 604 if (10 == fNumWarnings) { 605 ed << "\n ### G4SeltzerBergerModel Warnings stopped"; 606 } 607 G4Exception("G4SeltzerBergerModel::SampleScattering","em0044", 608 JustWarning, ed,""); 609 } 610 if (v >= vmax*rndm[1]) { 611 break; 612 } 613 } 614 return gammaEnergy; 615 } 616