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 // Author: Luciano Pandola 28 // 29 // History: 30 // -------- 31 // 23 Nov 2010 L Pandola First complete i 32 // 24 May 2011 L. Pandola Renamed (make de 33 // 13 Mar 2012 L. Pandola Updated the inte 34 // 18 Jul 2012 L. Pandola Migrate to the n 35 // now provides the 36 // 02 Oct 2013 L. Pandola Migrated to MT 37 // 17 Oct 2013 L. Pandola Partially revert 38 // kept as thread- 39 // 40 41 #include "G4PenelopeBremsstrahlungModel.hh" 42 #include "G4PhysicalConstants.hh" 43 #include "G4SystemOfUnits.hh" 44 #include "G4PenelopeBremsstrahlungFS.hh" 45 #include "G4PenelopeBremsstrahlungAngular.hh" 46 #include "G4ParticleDefinition.hh" 47 #include "G4MaterialCutsCouple.hh" 48 #include "G4ProductionCutsTable.hh" 49 #include "G4DynamicParticle.hh" 50 #include "G4Gamma.hh" 51 #include "G4Electron.hh" 52 #include "G4Positron.hh" 53 #include "G4PenelopeOscillatorManager.hh" 54 #include "G4PenelopeCrossSection.hh" 55 #include "G4PhysicsFreeVector.hh" 56 #include "G4PhysicsLogVector.hh" 57 #include "G4PhysicsTable.hh" 58 #include "G4Exp.hh" 59 60 //....oooOO0OOooo........oooOO0OOooo........oo 61 namespace { G4Mutex PenelopeBremsstrahlungMod 62 63 G4PenelopeBremsstrahlungModel::G4PenelopeBrems 64 const G4String& nam) 65 :G4VEmModel(nam),fParticleChange(nullptr),fP 66 fPenelopeFSHelper(nullptr),fPenelopeAngular 67 fXSTableElectron(nullptr),fXSTablePositron( 68 fIsInitialised(false),fLocalTable(false) 69 70 { 71 fIntrinsicLowEnergyLimit = 100.0*eV; 72 fIntrinsicHighEnergyLimit = 100.0*GeV; 73 nBins = 200; 74 75 if (part) 76 SetParticle(part); 77 78 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 79 // 80 fOscManager = G4PenelopeOscillatorManager::G 81 // 82 fVerboseLevel= 0; 83 // Verbosity scale: 84 // 0 = nothing 85 // 1 = warning for energy non-conservation 86 // 2 = details of energy budget 87 // 3 = calculation of cross sections, file o 88 // 4 = entering in methods 89 90 // Atomic deexcitation model activated by de 91 SetDeexcitationFlag(true); 92 93 } 94 95 //....oooOO0OOooo........oooOO0OOooo........oo 96 97 G4PenelopeBremsstrahlungModel::~G4PenelopeBrem 98 { 99 if (IsMaster() || fLocalTable) 100 { 101 ClearTables(); 102 if (fPenelopeFSHelper) 103 delete fPenelopeFSHelper; 104 } 105 // This is thread-local at the moment 106 if (fPenelopeAngular) 107 delete fPenelopeAngular; 108 } 109 110 //....oooOO0OOooo........oooOO0OOooo........oo 111 112 void G4PenelopeBremsstrahlungModel::Initialise 113 c 114 { 115 if (fVerboseLevel > 3) 116 G4cout << "Calling G4PenelopeBremsstrahlun 117 118 SetParticle(part); 119 120 if (IsMaster() && part == fParticle) 121 { 122 if (!fPenelopeFSHelper) 123 fPenelopeFSHelper = new G4PenelopeBremsstrah 124 if (!fPenelopeAngular) 125 fPenelopeAngular = new G4PenelopeBremsstrahl 126 //Clear and re-build the tables 127 ClearTables(); 128 129 //forces the cleaning of tables, in this 130 if (fPenelopeAngular) 131 fPenelopeAngular->Initialize(); 132 133 //Set the number of bins for the tables. 134 nBins = (std::size_t) (20*std::log10(Hig 135 nBins = std::max(nBins,(std::size_t)100) 136 fEnergyGrid = new G4PhysicsLogVector(Low 137 HighEner 138 nBins-1) 139 140 fXSTableElectron = new 141 std::map< std::pair<const G4Material*,G4doub 142 fXSTablePositron = new 143 std::map< std::pair<const G4Material*,G4doub 144 145 G4ProductionCutsTable* theCoupleTable = 146 G4ProductionCutsTable::GetProductionCutsTabl 147 148 //Build tables for all materials 149 for (G4int i=0;i<(G4int)theCoupleTable-> 150 { 151 const G4Material* theMat = 152 theCoupleTable->GetMaterialCutsCouple(i) 153 //Forces the building of the helper tables 154 fPenelopeFSHelper->BuildScaledXSTable(theM 155 fPenelopeAngular->PrepareTables(theMat,IsM 156 BuildXSTable(theMat,theCuts.at(i)); 157 158 } 159 160 if (fVerboseLevel > 2) { 161 G4cout << "Penelope Bremsstrahlung model v20 162 << "Energy range: " 163 << LowEnergyLimit() / keV << " keV - 164 << HighEnergyLimit() / GeV << " GeV." 165 << G4endl; 166 } 167 } 168 169 if(fIsInitialised) return; 170 fParticleChange = GetParticleChangeForLoss() 171 fIsInitialised = true; 172 } 173 174 //....oooOO0OOooo........oooOO0OOooo........oo 175 176 void G4PenelopeBremsstrahlungModel::Initialise 177 G4VEmModel *masterModel) 178 { 179 if (fVerboseLevel > 3) 180 G4cout << "Calling G4PenelopeBremsstrahlu 181 // 182 //Check that particle matches: one might hav 183 //for e+ and e-). 184 // 185 if (part == fParticle) 186 { 187 //Get the const table pointers from the 188 const G4PenelopeBremsstrahlungModel* the 189 static_cast<G4PenelopeBremsstrahlungModel*> 190 191 //Copy pointers to the data tables 192 fEnergyGrid = theModel->fEnergyGrid; 193 fXSTableElectron = theModel->fXSTableEle 194 fXSTablePositron = theModel->fXSTablePos 195 fPenelopeFSHelper = theModel->fPenelopeF 196 197 //created in each thread and initialized 198 if (!fPenelopeAngular) 199 fPenelopeAngular = new G4PenelopeBremsstrahl 200 //forces the cleaning of tables, in this 201 if (fPenelopeAngular) 202 fPenelopeAngular->Initialize(); 203 204 G4ProductionCutsTable* theCoupleTable = 205 G4ProductionCutsTable::GetProductionCutsTabl 206 //Build tables for all materials 207 for (G4int i=0;i<(G4int)theCoupleTable-> 208 { 209 const G4Material* theMat = 210 theCoupleTable->GetMaterialCutsCouple(i) 211 fPenelopeAngular->PrepareTables(theMat,IsM 212 } 213 214 //copy the data 215 nBins = theModel->nBins; 216 217 //Same verbosity for all workers, as the 218 fVerboseLevel = theModel->fVerboseLevel; 219 } 220 return; 221 } 222 223 //....oooOO0OOooo........oooOO0OOooo........oo 224 225 G4double G4PenelopeBremsstrahlungModel::CrossS 226 con 227 G4d 228 G4d 229 G4d 230 { 231 // 232 if (fVerboseLevel > 3) 233 G4cout << "Calling CrossSectionPerVolume() 234 235 SetupForMaterial(theParticle, material, ener 236 G4double crossPerMolecule = 0.; 237 238 G4PenelopeCrossSection* theXS = GetCrossSect 239 240 if (theXS) 241 crossPerMolecule = theXS->GetHardCrossSect 242 243 G4double atomDensity = material->GetTotNbOfA 244 G4double atPerMol = fOscManager->GetAtomsPe 245 246 if (fVerboseLevel > 3) 247 G4cout << "Material " << material->GetName 248 "atoms per molecule" << G4endl; 249 250 G4double moleculeDensity = 0.; 251 if (atPerMol) 252 moleculeDensity = atomDensity/atPerMol; 253 254 G4double crossPerVolume = crossPerMolecule*m 255 256 if (fVerboseLevel > 2) 257 { 258 G4cout << "G4PenelopeBremsstrahlungModel " 259 G4cout << "Mean free path for gamma emissi 260 energy/keV << " keV = " << (1./crossPerV 261 } 262 263 return crossPerVolume; 264 } 265 266 //....oooOO0OOooo........oooOO0OOooo........oo 267 268 //This is a dummy method. Never inkoved by the 269 //a warning if one tries to get Cross Sections 270 //G4EmCalculator. 271 G4double G4PenelopeBremsstrahlungModel::Comput 272 G4double, 273 G4double, 274 G4double, 275 G4double, 276 G4double) 277 { 278 G4cout << "*** G4PenelopeBremsstrahlungModel 279 G4cout << "Penelope Bremsstrahlung model v20 280 G4cout << "so the result is always zero. For 281 G4cout << "GetCrossSectionPerVolume() or Get 282 return 0; 283 } 284 285 //....oooOO0OOooo........oooOO0OOooo........oo 286 287 G4double G4PenelopeBremsstrahlungModel::Comput 288 const G4ParticleDefinition* 289 G4double kineticEnergy, 290 G4double cutEnergy) 291 { 292 if (fVerboseLevel > 3) 293 G4cout << "Calling ComputeDEDX() of G4Pene 294 295 G4PenelopeCrossSection* theXS = GetCrossSect 296 297 G4double sPowerPerMolecule = 0.0; 298 if (theXS) 299 sPowerPerMolecule = theXS->GetSoftStopping 300 301 G4double atomDensity = material->GetTotNbOfA 302 G4double atPerMol = fOscManager->GetAtomsPe 303 304 G4double moleculeDensity = 0.; 305 if (atPerMol) 306 moleculeDensity = atomDensity/atPerMol; 307 308 G4double sPowerPerVolume = sPowerPerMolecule 309 310 if (fVerboseLevel > 2) 311 { 312 G4cout << "G4PenelopeBremsstrahlungModel 313 G4cout << "Stopping power < " << cutEner 314 kineticEnergy/keV << " keV = " << 315 sPowerPerVolume/(keV/mm) << " keV/mm" 316 } 317 return sPowerPerVolume; 318 } 319 320 //....oooOO0OOooo........oooOO0OOooo........oo 321 322 void G4PenelopeBremsstrahlungModel::SampleSeco 323 const G4MaterialCutsCouple* 324 const G4DynamicParticle* aDy 325 G4double cutG, 326 G4double) 327 { 328 if (fVerboseLevel > 3) 329 G4cout << "Calling SampleSecondaries() of 330 331 G4double kineticEnergy = aDynamicParticle->G 332 const G4Material* material = couple->GetMate 333 334 if (kineticEnergy <= fIntrinsicLowEnergyLimi 335 { 336 fParticleChange->SetProposedKineticEnergy 337 fParticleChange->ProposeLocalEnergyDeposi 338 return ; 339 } 340 341 G4ParticleMomentum particleDirection0 = aDyn 342 //This is the momentum 343 G4ThreeVector initialMomentum = aDynamicPar 344 345 //Not enough energy to produce a secondary! 346 if (kineticEnergy < cutG) 347 return; 348 349 if (fVerboseLevel > 3) 350 G4cout << "Going to sample gamma energy fo 351 "energy = " << kineticEnergy/keV << ", c 352 353 //Sample gamma's energy according to the sp 354 G4double gammaEnergy = 355 fPenelopeFSHelper->SampleGammaEnergy(kinet 356 357 if (fVerboseLevel > 3) 358 G4cout << "Sampled gamma energy: " << gamm 359 360 //Now sample the direction for the Gamma. No 361 //Z is unused here, I plug 0. The informatio 362 G4ThreeVector gammaDirection1 = 363 fPenelopeAngular->SampleDirection(aDynami 364 365 if (fVerboseLevel > 3) 366 G4cout << "Sampled cosTheta for e-: " << g 367 368 G4double residualPrimaryEnergy = kineticEner 369 if (residualPrimaryEnergy < 0) 370 { 371 //Ok we have a problem, all energy goes 372 gammaEnergy += residualPrimaryEnergy; 373 residualPrimaryEnergy = 0.0; 374 } 375 376 //Produce final state according to momentum 377 G4ThreeVector particleDirection1 = initialMo 378 particleDirection1 = particleDirection1.unit 379 380 //Update the primary particle 381 if (residualPrimaryEnergy > 0.) 382 { 383 fParticleChange->ProposeMomentumDirectio 384 fParticleChange->SetProposedKineticEnerg 385 } 386 else 387 { 388 fParticleChange->SetProposedKineticEnerg 389 } 390 391 //Now produce the photon 392 G4DynamicParticle* theGamma = new G4DynamicP 393 394 395 fvect->push_back(theGamma); 396 397 if (fVerboseLevel > 1) 398 { 399 G4cout << "----------------------------- 400 G4cout << "Energy balance from G4Penelop 401 G4cout << "Incoming primary energy: " << 402 G4cout << "----------------------------- 403 G4cout << "Outgoing primary energy: " << 404 G4cout << "Bremsstrahlung photon " << ga 405 G4cout << "Total final state: " << (resi 406 << " keV" << G4endl; 407 G4cout << "----------------------------- 408 } 409 410 if (fVerboseLevel > 0) 411 { 412 G4double energyDiff = std::fabs(residual 413 if (energyDiff > 0.05*keV) 414 G4cout << "Warning from G4PenelopeBrem 415 << 416 (residualPrimaryEnergy+gammaEnergy)/ 417 " keV (final) vs. " << 418 kineticEnergy/keV << " keV (initial) 419 } 420 return; 421 } 422 423 //....oooOO0OOooo........oooOO0OOooo........oo 424 425 void G4PenelopeBremsstrahlungModel::ClearTable 426 { 427 if (!IsMaster() && !fLocalTable) 428 //Should not be here! 429 G4Exception("G4PenelopeBremsstrahlungModel 430 "em0100",FatalException,"Worker thread in 431 432 if (fXSTableElectron) 433 { 434 for (auto& item : (*fXSTableElectron)) 435 delete item.second; 436 delete fXSTableElectron; 437 fXSTableElectron = nullptr; 438 } 439 if (fXSTablePositron) 440 { 441 for (auto& item : (*fXSTablePositron)) 442 delete item.second; 443 delete fXSTablePositron; 444 fXSTablePositron = nullptr; 445 } 446 /* 447 if (fEnergyGrid) 448 delete fEnergyGrid; 449 */ 450 if (fPenelopeFSHelper) 451 fPenelopeFSHelper->ClearTables(IsMaster()) 452 453 if (fVerboseLevel > 2) 454 G4cout << "G4PenelopeBremsstrahlungModel: 455 456 return; 457 } 458 459 //....oooOO0OOooo........oooOO0OOooo........oo 460 461 G4double G4PenelopeBremsstrahlungModel::MinEne 462 const G4MaterialCutsCouple* 463 { 464 return fIntrinsicLowEnergyLimit; 465 } 466 467 //....oooOO0OOooo........oooOO0OOooo........oo 468 469 void G4PenelopeBremsstrahlungModel::BuildXSTab 470 { 471 if (!IsMaster() && !fLocalTable) 472 //Should not be here! 473 G4Exception("G4PenelopeBremsstrahlungModel 474 "em0100",FatalException,"Worker thread in 475 476 //The key of the map 477 std::pair<const G4Material*,G4double> theKey 478 479 //The table already exists 480 if (fXSTableElectron->count(theKey) && fXSTa 481 return; 482 483 // 484 //This method fills the G4PenelopeCrossSecti 485 //and for the given material/cut couple. 486 //Equivalent of subroutines EBRaT and PINaT 487 // 488 if (fVerboseLevel > 2) 489 { 490 G4cout << "G4PenelopeBremsstrahlungModel 491 G4cout << "for e+/e- in " << mat->GetNam 492 cut/keV << " keV " << G4endl; 493 } 494 495 //Tables have been already created (checked 496 if (fEnergyGrid->GetVectorLength() != nBins) 497 { 498 G4ExceptionDescription ed; 499 ed << "Energy Grid looks not initialized 500 ed << nBins << " " << fEnergyGrid->GetVe 501 G4Exception("G4PenelopeBremsstrahlungMod 502 "em2016",FatalException,ed); 503 } 504 505 G4PenelopeCrossSection* XSEntry = new G4Pene 506 G4PenelopeCrossSection* XSEntry2 = new G4Pen 507 const G4PhysicsTable* table = fPenelopeFSHel 508 509 //loop on the energy grid 510 for (std::size_t bin=0;bin<nBins;++bin) 511 { 512 G4double energy = fEnergyGrid->GetLowEd 513 G4double XH0=0, XH1=0, XH2=0; 514 G4double XS0=0, XS1=0, XS2=0; 515 516 //Global xs factor 517 G4double fact = fPenelopeFSHelper->GetE 518 ((energy+electron_mass_c2)*(energy+electron 519 (energy*(energy+2.0*electron_mass_c2))); 520 521 G4double restrictedCut = cut/energy; 522 523 //Now I need the dSigma/dX profile - in 524 //the 32-point x grid. Interpolation is 525 std::size_t nBinsX = fPenelopeFSHelper- 526 G4double* tempData = new G4double[nBins 527 G4double logene = G4Log(energy); 528 for (std::size_t ix=0;ix<nBinsX;++ix) 529 { 530 //find dSigma/dx for the given E. X belon 531 G4double val = (*table)[ix]->Value(logene 532 tempData[ix] = G4Exp(val); //back to the 533 } 534 535 G4double XH0A = 0.; 536 if (restrictedCut <= 1) //calculate onl 537 XH0A = fPenelopeFSHelper->GetMomentumIntegr 538 fPenelopeFSHelper->GetMomentumIntegral(te 539 G4double XS1A = fPenelopeFSHelper->GetM 540 restrictedCut,0); 541 G4double XS2A = fPenelopeFSHelper->GetM 542 restrictedCut,1); 543 G4double XH1A=0, XH2A=0; 544 if (restrictedCut <=1) 545 { 546 XH1A = fPenelopeFSHelper->GetMomentumInte 547 XS1A; 548 XH2A = fPenelopeFSHelper->GetMomentumInte 549 XS2A; 550 } 551 delete[] tempData; 552 553 XH0 = XH0A*fact; 554 XS1 = XS1A*fact*energy; 555 XH1 = XH1A*fact*energy; 556 XS2 = XS2A*fact*energy*energy; 557 XH2 = XH2A*fact*energy*energy; 558 559 XSEntry->AddCrossSectionPoint(bin,energ 560 561 //take care of positrons 562 G4double posCorrection = GetPositronXSC 563 XSEntry2->AddCrossSectionPoint(bin,ener 564 XH1*posCorrection, 565 XH2*posCorrection, 566 XS0, 567 XS1*posCorrection, 568 XS2*posCorrection); 569 } 570 571 //Insert in the appropriate table 572 fXSTableElectron->insert(std::make_pair(theK 573 fXSTablePositron->insert(std::make_pair(theK 574 575 return; 576 } 577 578 579 //....oooOO0OOooo........oooOO0OOooo........oo 580 581 G4PenelopeCrossSection* 582 G4PenelopeBremsstrahlungModel::GetCrossSection 583 const G4Material* mat, 584 G4double cut) 585 { 586 if (part != G4Electron::Electron() && part ! 587 { 588 G4ExceptionDescription ed; 589 ed << "Invalid particle: " << part->GetP 590 G4Exception("G4PenelopeBremsstrahlungMod 591 "em0001",FatalException,ed); 592 return nullptr; 593 } 594 595 if (part == G4Electron::Electron()) 596 { 597 //Either Initialize() was not called, or 598 //not invoked 599 if (!fXSTableElectron) 600 { 601 //create a **thread-local** version of the 602 //Unit Tests 603 G4String excep = "The Cross Section 604 G4Exception("G4PenelopeBremsstrahlun 605 "em2013",JustWarning,excep); 606 fLocalTable = true; 607 fXSTableElectron = new 608 std::map< std::pair<const G4Material*,G4 609 if (!fEnergyGrid) 610 fEnergyGrid = new G4PhysicsLogVector(Low 611 HighEnergyLimit(), 612 nBins-1); //one hidden bin is adde 613 if (!fPenelopeFSHelper) 614 fPenelopeFSHelper = new G4PenelopeBremss 615 } 616 std::pair<const G4Material*,G4double> th 617 if (fXSTableElectron->count(theKey)) //t 618 return fXSTableElectron->find(theKey)- 619 else 620 { 621 //If we are here, it means that Initialize 622 //not filled up. This can happen in a Unit 623 if (fVerboseLevel > 0) 624 { 625 //G4Exception (warning) is issued only 626 G4ExceptionDescription ed; 627 ed << "Unable to find e- table for " < 628 << cut/keV << " keV " << G4endl; 629 ed << "This can happen only in Unit Te 630 G4Exception("G4PenelopeBremsstrahlungM 631 "em2009",JustWarning,ed); 632 } 633 //protect file reading via autolock 634 G4AutoLock lock(&PenelopeBremsstrahlungMod 635 fPenelopeFSHelper->BuildScaledXSTabl 636 BuildXSTable(mat,cut); 637 lock.unlock(); 638 //now it should be ok 639 return fXSTableElectron->find(theKey)->sec 640 } 641 } 642 if (part == G4Positron::Positron()) 643 { 644 //Either Initialize() was not called, or 645 //not invoked 646 if (!fXSTablePositron) 647 { 648 G4String excep = "The Cross Section Table 649 G4Exception("G4PenelopeBremsstrahlun 650 "em2013",JustWarning,excep); 651 fLocalTable = true; 652 fXSTablePositron = new 653 std::map< std::pair<const G4Material*,G4 654 if (!fEnergyGrid) 655 fEnergyGrid = new G4PhysicsLogVector(Low 656 HighEnergyLimit(), 657 nBins-1); //one hidden bin is adde 658 if (!fPenelopeFSHelper) 659 fPenelopeFSHelper = new G4Penelope 660 } 661 std::pair<const G4Material*,G4double> th 662 if (fXSTablePositron->count(theKey)) //t 663 return fXSTablePositron->find(theKey)- 664 else 665 { 666 //If we are here, it means that Initialize 667 //not filled up. This can happen in a Unit 668 if (fVerboseLevel > 0) 669 { 670 //Issue a G4Exception (warning) only i 671 G4ExceptionDescription ed; 672 ed << "Unable to find e+ table for " < 673 << cut/keV << " keV " << G4endl; 674 ed << "This can happen only in Unit Te 675 G4Exception("G4PenelopeBremsstrahlungM 676 "em2009",JustWarning,ed); 677 } 678 //protect file reading via autolock 679 G4AutoLock lock(&PenelopeBremsstrahlungMod 680 fPenelopeFSHelper->BuildScaledXSTabl 681 BuildXSTable(mat,cut); 682 lock.unlock(); 683 //now it should be ok 684 return fXSTablePositron->find(theKey)->sec 685 } 686 } 687 return nullptr; 688 } 689 690 //....oooOO0OOooo........oooOO0OOooo........oo 691 692 G4double G4PenelopeBremsstrahlungModel::GetPos 693 G4double energy) 694 { 695 //The electron-to-positron correction factor 696 //radiative stopping powers for positrons an 697 //by Kim et al. (1986) (cf. Berger and Seltz 698 //analytical approximation which reproduces 699 //accuracy 700 G4double t=G4Log(1.0+1e6*energy/ 701 (electron_mass_c2*fPenelopeFSHelper- 702 G4double corr = 1.0-G4Exp(-t*(1.2359e-1-t*(6 703 (3.1516e-2-t*(7.7446e-3-t*(1.0595 704 (7.0568e-5-t* 705 1.8080e-6))))))); 706 return corr; 707 } 708 709 //....oooOO0OOooo........oooOO0OOooo........oo 710 711 void G4PenelopeBremsstrahlungModel::SetParticl 712 { 713 if(!fParticle) { 714 fParticle = p; 715 } 716 } 717