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 // Author: Luciano Pandola 28 // 29 // History: 30 // -------- 31 // 23 Nov 2010 L Pandola First complete implementation 32 // 02 May 2011 L.Pandola Remove dependency on CLHEP::HepMatrix 33 // 24 May 2011 L.Pandola Renamed (make v2008 as default Penelope) 34 // 03 Oct 2013 L.Pandola Migration to MT 35 // 30 Oct 2013 L.Pandola Use G4Cache to avoid new/delete of the 36 // data vector on the fly in SampleGammaEnergy() 37 // 38 39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 40 #include "G4PenelopeBremsstrahlungFS.hh" 41 #include "G4PhysicsFreeVector.hh" 42 #include "G4PhysicsTable.hh" 43 #include "G4Material.hh" 44 #include "Randomize.hh" 45 #include "G4AutoDelete.hh" 46 #include "G4Exp.hh" 47 #include "G4PhysicalConstants.hh" 48 #include "G4SystemOfUnits.hh" 49 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 51 52 G4PenelopeBremsstrahlungFS::G4PenelopeBremsstrahlungFS(G4int verbosity) : 53 fReducedXSTable(nullptr),fEffectiveZSq(nullptr),fSamplingTable(nullptr), 54 fPBcut(nullptr),fVerbosity(verbosity) 55 { 56 fCache.Put(0); 57 G4double tempvector[fNBinsX] = 58 {1.0e-12,0.025e0,0.05e0,0.075e0,0.1e0,0.15e0,0.2e0,0.25e0, 59 0.3e0,0.35e0,0.4e0,0.45e0,0.5e0,0.55e0,0.6e0,0.65e0,0.7e0, 60 0.75e0,0.8e0,0.85e0,0.9e0,0.925e0,0.95e0,0.97e0,0.99e0, 61 0.995e0,0.999e0,0.9995e0,0.9999e0,0.99995e0,0.99999e0,1.0e0}; 62 63 for (std::size_t ix=0;ix<fNBinsX;++ix) 64 theXGrid[ix] = tempvector[ix]; 65 66 for (std::size_t i=0;i<fNBinsE;++i) 67 theEGrid[i] = 0.; 68 69 fElementData = new std::map<G4int,G4DataVector*>; 70 } 71 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 73 74 G4PenelopeBremsstrahlungFS::~G4PenelopeBremsstrahlungFS() 75 { 76 ClearTables(); 77 78 //The G4Physics*Vector pointers contained in the fCache are automatically deleted by 79 //the G4AutoDelete so there is no need to take care of them manually 80 81 //Clear manually fElementData 82 if (fElementData) 83 { 84 for (auto& item : (*fElementData)) 85 delete item.second; 86 delete fElementData; 87 fElementData = nullptr; 88 } 89 } 90 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo... 92 93 94 void G4PenelopeBremsstrahlungFS::ClearTables(G4bool isMaster) 95 { 96 //Just to check 97 if (!isMaster) 98 G4Exception("G4PenelopeBremsstrahlungFS::ClearTables()", 99 "em0100",FatalException,"Worker thread in this method"); 100 101 if (fReducedXSTable) 102 { 103 for (auto& item : (*fReducedXSTable)) 104 { 105 G4PhysicsTable* tab = item.second; 106 tab->clearAndDestroy(); 107 delete tab; 108 } 109 fReducedXSTable->clear(); 110 delete fReducedXSTable; 111 fReducedXSTable = nullptr; 112 } 113 114 if (fSamplingTable) 115 { 116 for (auto& item : (*fSamplingTable)) 117 { 118 G4PhysicsTable* tab = item.second; 119 tab->clearAndDestroy(); 120 delete tab; 121 } 122 fSamplingTable->clear(); 123 delete fSamplingTable; 124 fSamplingTable = nullptr; 125 } 126 if (fPBcut) 127 { 128 /* 129 std::map< std::pair<const G4Material*,G4double> ,G4PhysicsFreeVector*>::iterator kk; 130 for (kk=fPBcut->begin(); kk != fPBcut->end(); kk++) 131 delete kk->second; 132 */ 133 delete fPBcut; 134 fPBcut = nullptr; 135 } 136 137 if (fEffectiveZSq) 138 { 139 delete fEffectiveZSq; 140 fEffectiveZSq = nullptr; 141 } 142 143 return; 144 } 145 146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 147 148 G4double G4PenelopeBremsstrahlungFS::GetEffectiveZSquared(const G4Material* material) const 149 { 150 if (!fEffectiveZSq) 151 { 152 G4ExceptionDescription ed; 153 ed << "The container for the <Z^2> values is not initialized" << G4endl; 154 G4Exception("G4PenelopeBremsstrahlungFS::GetEffectiveZSquared()", 155 "em2007",FatalException,ed); 156 return 0; 157 } 158 //found in the table: return it 159 if (fEffectiveZSq->count(material)) 160 return fEffectiveZSq->find(material)->second; 161 else 162 { 163 G4ExceptionDescription ed; 164 ed << "The value of <Z^2> is not properly set for material " << 165 material->GetName() << G4endl; 166 //requires running of BuildScaledXSTable() 167 G4Exception("G4PenelopeBremsstrahlungFS::GetEffectiveZSquared()", 168 "em2008",FatalException,ed); 169 } 170 return 0; 171 } 172 173 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 174 175 void G4PenelopeBremsstrahlungFS::BuildScaledXSTable(const G4Material* material, 176 G4double cut,G4bool isMaster) 177 { 178 //Corresponds to subroutines EBRaW and EBRaR of PENELOPE 179 /* 180 This method generates the table of the scaled energy-loss cross section from 181 bremsstrahlung emission for the given material. Original data are read from 182 file. The table is normalized according to the Berger-Seltzer cross section. 183 */ 184 185 //Just to check 186 if (!isMaster) 187 G4Exception("G4PenelopeBremsstrahlungFS::BuildScaledXSTable()", 188 "em0100",FatalException,"Worker thread in this method"); 189 190 if (fVerbosity > 2) 191 { 192 G4cout << "Entering in G4PenelopeBremsstrahlungFS::BuildScaledXSTable for " << 193 material->GetName() << G4endl; 194 G4cout << "Threshold = " << cut/keV << " keV, isMaster= " << isMaster << 195 G4endl; 196 } 197 198 //This method should be accessed by the master only 199 if (!fSamplingTable) 200 fSamplingTable = 201 new std::map< std::pair<const G4Material*,G4double> , G4PhysicsTable*>; 202 if (!fPBcut) 203 fPBcut = 204 new std::map< std::pair<const G4Material*,G4double> , G4PhysicsFreeVector* >; 205 206 //check if the container exists (if not, create it) 207 if (!fReducedXSTable) 208 fReducedXSTable = new std::map< std::pair<const G4Material*,G4double> , 209 G4PhysicsTable*>; 210 if (!fEffectiveZSq) 211 fEffectiveZSq = new std::map<const G4Material*,G4double>; 212 213 //********************************************************************* 214 //Determine the equivalent atomic number <Z^2> 215 //********************************************************************* 216 std::vector<G4double> *StechiometricFactors = new std::vector<G4double>; 217 std::size_t nElements = material->GetNumberOfElements(); 218 const G4ElementVector* elementVector = material->GetElementVector(); 219 const G4double* fractionVector = material->GetFractionVector(); 220 for (std::size_t i=0;i<nElements;i++) 221 { 222 G4double fraction = fractionVector[i]; 223 G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole); 224 StechiometricFactors->push_back(fraction/atomicWeigth); 225 } 226 //Find max 227 G4double MaxStechiometricFactor = 0.; 228 for (std::size_t i=0;i<nElements;i++) 229 { 230 if ((*StechiometricFactors)[i] > MaxStechiometricFactor) 231 MaxStechiometricFactor = (*StechiometricFactors)[i]; 232 } 233 //Normalize 234 for (std::size_t i=0;i<nElements;i++) 235 (*StechiometricFactors)[i] /= MaxStechiometricFactor; 236 237 G4double sumz2 = 0; 238 G4double sums = 0; 239 for (std::size_t i=0;i<nElements;i++) 240 { 241 G4double Z = (*elementVector)[i]->GetZ(); 242 sumz2 += (*StechiometricFactors)[i]*Z*Z; 243 sums += (*StechiometricFactors)[i]; 244 } 245 G4double ZBR2 = sumz2/sums; 246 247 fEffectiveZSq->insert(std::make_pair(material,ZBR2)); 248 249 //********************************************************************* 250 // loop on elements and read data files 251 //********************************************************************* 252 G4DataVector* tempData = new G4DataVector(fNBinsE); 253 G4DataVector* tempMatrix = new G4DataVector(fNBinsE*fNBinsX,0.); 254 255 for (std::size_t iel=0;iel<nElements;iel++) 256 { 257 G4double Z = (*elementVector)[iel]->GetZ(); 258 G4int iZ = (G4int) Z; 259 G4double wgt = (*StechiometricFactors)[iel]*Z*Z/ZBR2; 260 261 //the element is not already loaded 262 if (!fElementData->count(iZ)) 263 { 264 ReadDataFile(iZ); 265 if (!fElementData->count(iZ)) 266 { 267 G4ExceptionDescription ed; 268 ed << "Error in G4PenelopeBremsstrahlungFS::BuildScaledXSTable" << G4endl; 269 ed << "Unable to retrieve data for element " << iZ << G4endl; 270 G4Exception("G4PenelopeBremsstrahlungFS::BuildScaledXSTable()", 271 "em2009",FatalException,ed); 272 } 273 } 274 275 G4DataVector* atomData = fElementData->find(iZ)->second; 276 277 for (std::size_t ie=0;ie<fNBinsE;++ie) 278 { 279 (*tempData)[ie] += wgt*(*atomData)[ie*(fNBinsX+1)+fNBinsX]; //last column contains total XS 280 for (std::size_t ix=0;ix<fNBinsX;++ix) 281 (*tempMatrix)[ie*fNBinsX+ix] += wgt*(*atomData)[ie*(fNBinsX+1)+ix]; 282 } 283 } 284 285 //********************************************************************* 286 // the total energy loss spectrum is re-normalized to reproduce the total 287 // scaled cross section of Berger and Seltzer 288 //********************************************************************* 289 for (std::size_t ie=0;ie<fNBinsE;++ie) 290 { 291 //for each energy, calculate integral of dSigma/dx over dx 292 G4double* tempData2 = new G4double[fNBinsX]; 293 for (std::size_t ix=0;ix<fNBinsX;++ix) 294 tempData2[ix] = (*tempMatrix)[ie*fNBinsX+ix]; 295 G4double rsum = GetMomentumIntegral(tempData2,1.0,0); 296 delete[] tempData2; 297 G4double fact = millibarn*(theEGrid[ie]+electron_mass_c2)*(1./fine_structure_const)/ 298 (classic_electr_radius*classic_electr_radius*(theEGrid[ie]+2.0*electron_mass_c2)); 299 G4double fnorm = (*tempData)[ie]/(rsum*fact); 300 G4double TST = 100.*std::fabs(fnorm-1.0); 301 if (TST > 1.0) 302 { 303 G4ExceptionDescription ed; 304 ed << "G4PenelopeBremsstrahlungFS. Corrupted data files?" << G4endl; 305 G4cout << "TST= " << TST << "; fnorm = " << fnorm << G4endl; 306 G4cout << "rsum = " << rsum << G4endl; 307 G4cout << "fact = " << fact << G4endl; 308 G4cout << ie << " " << theEGrid[ie]/keV << " " << (*tempData)[ie]/barn << G4endl; 309 G4Exception("G4PenelopeBremsstrahlungFS::BuildScaledXSTable()", 310 "em2010",FatalException,ed); 311 } 312 for (std::size_t ix=0;ix<fNBinsX;++ix) 313 (*tempMatrix)[ie*fNBinsX+ix] *= fnorm; 314 } 315 316 //********************************************************************* 317 // create and fill the tables 318 //********************************************************************* 319 G4PhysicsTable* thePhysicsTable = new G4PhysicsTable(); 320 // the table will contain 32 G4PhysicsFreeVectors with different 321 // values of x. Each of the G4PhysicsFreeVectors has a profile of 322 // log(XS) vs. log(E) 323 324 //reserve space of the vectors. Everything is log-log 325 //I add one extra "fake" point at low energy, since the Penelope 326 //table starts at 1 keV 327 for (std::size_t i=0;i<fNBinsX;++i) 328 thePhysicsTable->push_back(new G4PhysicsFreeVector(fNBinsE+1)); 329 330 for (std::size_t ix=0;ix<fNBinsX;++ix) 331 { 332 G4PhysicsFreeVector* theVec = 333 (G4PhysicsFreeVector*) ((*thePhysicsTable)[ix]); 334 for (std::size_t ie=0;ie<fNBinsE;++ie) 335 { 336 G4double logene = G4Log(theEGrid[ie]); 337 G4double aValue = (*tempMatrix)[ie*fNBinsX+ix]; 338 if (aValue < 1e-20*millibarn) //protection against log(0) 339 aValue = 1e-20*millibarn; 340 theVec->PutValues(ie+1,logene,G4Log(aValue)); 341 } 342 //Add fake point at 1 eV using an extrapolation with the derivative 343 //at the first valid point (Penelope approach) 344 G4double derivative = ((*theVec)[2]-(*theVec)[1])/(theVec->Energy(2) - theVec->Energy(1)); 345 G4double log1eV = G4Log(1*eV); 346 G4double val1eV = (*theVec)[1]+derivative*(log1eV-theVec->Energy(1)); 347 //fake point at very low energy 348 theVec->PutValues(0,log1eV,val1eV); 349 } 350 std::pair<const G4Material*,G4double> theKey = std::make_pair(material,cut); 351 fReducedXSTable->insert(std::make_pair(theKey,thePhysicsTable)); 352 353 delete StechiometricFactors; 354 delete tempData; 355 delete tempMatrix; 356 357 //Do here also the initialization of the energy sampling 358 if (!(fSamplingTable->count(theKey))) 359 InitializeEnergySampling(material,cut); 360 361 return; 362 } 363 364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 365 366 void G4PenelopeBremsstrahlungFS::ReadDataFile(G4int Z) 367 { 368 const char* path = G4FindDataDir("G4LEDATA"); 369 if (!path) 370 { 371 G4String excep = "G4PenelopeBremsstrahlungFS - G4LEDATA environment variable not set!"; 372 G4Exception("G4PenelopeBremsstrahlungFS::ReadDataFile()", 373 "em0006",FatalException,excep); 374 return; 375 } 376 /* 377 Read the cross section file 378 */ 379 std::ostringstream ost; 380 if (Z>9) 381 ost << path << "/penelope/bremsstrahlung/pdebr" << Z << ".p08"; 382 else 383 ost << path << "/penelope/bremsstrahlung/pdebr0" << Z << ".p08"; 384 std::ifstream file(ost.str().c_str()); 385 if (!file.is_open()) 386 { 387 G4String excep = "G4PenelopeBremsstrahlungFS - data file " + 388 G4String(ost.str()) + " not found!"; 389 G4Exception("G4PenelopeBremsstrahlungFS::ReadDataFile()", 390 "em0003",FatalException,excep); 391 return; 392 } 393 394 G4int readZ =0; 395 file >> readZ; 396 397 //check the right file is opened. 398 if (readZ != Z) 399 { 400 G4ExceptionDescription ed; 401 ed << "Corrupted data file for Z=" << Z << G4endl; 402 G4Exception("G4PenelopeBremsstrahlungFS::ReadDataFile()", 403 "em0005",FatalException,ed); 404 return; 405 } 406 407 G4DataVector* theMatrix = new G4DataVector(fNBinsE*(fNBinsX+1),0.); //initialized with zeros 408 409 for (std::size_t ie=0;ie<fNBinsE;++ie) 410 { 411 G4double myDouble = 0; 412 file >> myDouble; //energy (eV) 413 if (!theEGrid[ie]) //fill only the first time 414 theEGrid[ie] = myDouble*eV; 415 // 416 for (std::size_t ix=0;ix<fNBinsX;++ix) 417 { 418 file >> myDouble; 419 (*theMatrix)[ie*(fNBinsX+1)+ix] = myDouble*millibarn; 420 } 421 file >> myDouble; //total cross section 422 (*theMatrix)[ie*(fNBinsX+1)+fNBinsX] = myDouble*millibarn; 423 } 424 425 if (fElementData) 426 fElementData->insert(std::make_pair(Z,theMatrix)); 427 else 428 delete theMatrix; 429 file.close(); 430 return; 431 } 432 433 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 434 435 G4double G4PenelopeBremsstrahlungFS::GetMomentumIntegral(G4double* y, 436 G4double xup,G4int momOrder) const 437 //x is always the gridX 438 { 439 //Corresponds to the function RLMOM of Penelope 440 //This method performs the calculation of the integral of (x^momOrder)*y over the interval 441 //from x[0] to xup, obtained by linear interpolation on a table of y. 442 //The independent variable is assumed to take positive values only. 443 // 444 std::size_t size = fNBinsX; 445 const G4double eps = 1e-35; 446 447 //Check that the call is valid 448 if (momOrder<-1 || size<2 || theXGrid[0]<0) 449 { 450 G4Exception("G4PenelopeBremsstrahlungFS::GetMomentumIntegral()", 451 "em2011",FatalException,"Invalid call"); 452 } 453 454 for (std::size_t i=1;i<size;++i) 455 { 456 if (theXGrid[i]<0 || theXGrid[i]<theXGrid[i-1]) 457 { 458 G4ExceptionDescription ed; 459 ed << "Invalid call for bin " << i << G4endl; 460 G4Exception("G4PenelopeBremsstrahlungFS::GetMomentumIntegral()", 461 "em2012",FatalException,ed); 462 } 463 } 464 465 //Compute the integral 466 G4double result = 0; 467 if (xup < theXGrid[0]) 468 return result; 469 G4bool loopAgain = true; 470 G4double xt = std::min(xup,theXGrid[size-1]); 471 G4double xtc = 0; 472 for (std::size_t i=0;i<size-1;++i) 473 { 474 G4double x1 = std::max(theXGrid[i],eps); 475 G4double y1 = y[i]; 476 G4double x2 = std::max(theXGrid[i+1],eps); 477 G4double y2 = y[i+1]; 478 if (xt < x2) 479 { 480 xtc = xt; 481 loopAgain = false; 482 } 483 else 484 xtc = x2; 485 G4double dx = x2-x1; 486 G4double dy = y2-y1; 487 G4double ds = 0; 488 if (std::fabs(dx)>1e-14*std::fabs(dy)) 489 { 490 G4double b=dy/dx; 491 G4double a=y1-b*x1; 492 if (momOrder == -1) 493 ds = a*G4Log(xtc/x1)+b*(xtc-x1); 494 else if (momOrder == 0) //speed it up, not using pow() 495 ds = a*(xtc-x1) + 0.5*b*(xtc*xtc-x1*x1); 496 else 497 ds = a*(std::pow(xtc,momOrder+1)-std::pow(x1,momOrder+1))/((G4double) (momOrder + 1)) 498 + b*(std::pow(xtc,momOrder+2)-std::pow(x1,momOrder+2))/((G4double) (momOrder + 2)); 499 } 500 else 501 ds = 0.5*(y1+y2)*(xtc-x1)*std::pow(xtc,momOrder); 502 result += ds; 503 if (!loopAgain) 504 return result; 505 } 506 return result; 507 } 508 509 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 510 511 const G4PhysicsTable* G4PenelopeBremsstrahlungFS::GetScaledXSTable(const G4Material* mat, 512 const G4double cut) const 513 { 514 //check if it already contains the entry 515 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut); 516 517 if (!(fReducedXSTable->count(theKey))) 518 { 519 G4Exception("G4PenelopeBremsstrahlungFS::GetScaledXSTable()", 520 "em2013",FatalException,"Unable to retrieve the cross section table"); 521 } 522 523 return fReducedXSTable->find(theKey)->second; 524 } 525 526 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 527 528 void G4PenelopeBremsstrahlungFS::InitializeEnergySampling(const G4Material* material, 529 G4double cut) 530 { 531 if (fVerbosity > 2) 532 G4cout << "Entering in G4PenelopeBremsstrahlungFS::InitializeEnergySampling() for " << 533 material->GetName() << G4endl; 534 535 //This method should be accessed by the master only 536 std::pair<const G4Material*,G4double> theKey = std::make_pair(material,cut); 537 538 G4PhysicsTable* thePhysicsTable = new G4PhysicsTable(); 539 // the table will contain 57 G4PhysicsFreeVectors with different 540 // values of E. 541 G4PhysicsFreeVector* thePBvec = new G4PhysicsFreeVector(fNBinsE); 542 543 //I reserve space of the vectors. 544 for (std::size_t i=0;i<fNBinsE;++i) 545 thePhysicsTable->push_back(new G4PhysicsFreeVector(fNBinsX)); 546 547 //Retrieve the table. Must already exist at this point, because this 548 //method is invoked by GetScaledXSTable() 549 if (!(fReducedXSTable->count(theKey))) 550 G4Exception("G4PenelopeBremsstrahlungFS::InitializeEnergySampling()", 551 "em2013",FatalException,"Unable to retrieve the cross section table"); 552 G4PhysicsTable* theTableReduced = fReducedXSTable->find(theKey)->second; 553 554 for (std::size_t ie=0;ie<fNBinsE;++ie) 555 { 556 G4PhysicsFreeVector* theVec = 557 (G4PhysicsFreeVector*) ((*thePhysicsTable)[ie]); 558 //Fill the table 559 G4double value = 0; //first value 560 theVec->PutValues(0,theXGrid[0],value); 561 for (std::size_t ix=1;ix<fNBinsX;++ix) 562 { 563 //Here calculate the cumulative distribution 564 // int_{0}^{x} dSigma(x',E)/dx' (1/x') dx' 565 G4PhysicsFreeVector* v1 = (G4PhysicsFreeVector*) (*theTableReduced)[ix-1]; 566 G4PhysicsFreeVector* v2 = (G4PhysicsFreeVector*) (*theTableReduced)[ix]; 567 568 G4double x1=std::max(theXGrid[ix-1],1.0e-35); 569 //Remember: the table fReducedXSTable has a fake first point in energy 570 //so, it contains one more bin than fNBinsE. 571 G4double y1=G4Exp((*v1)[ie+1]); 572 G4double x2=std::max(theXGrid[ix],1.0e-35); 573 G4double y2=G4Exp((*v2)[ie+1]); 574 G4double B = (y2-y1)/(x2-x1); 575 G4double A = y1-B*x1; 576 G4double dS = A*G4Log(x2/x1)+B*(x2-x1); 577 value += dS; 578 theVec->PutValues(ix,theXGrid[ix],value); 579 } 580 //fill the PB vector 581 G4double xc = cut/theEGrid[ie]; 582 //Fill a temp data vector 583 G4double* tempData = new G4double[fNBinsX]; 584 for (std::size_t ix=0;ix<fNBinsX;++ix) 585 { 586 G4PhysicsFreeVector* vv = (G4PhysicsFreeVector*) (*theTableReduced)[ix]; 587 tempData[ix] = G4Exp((*vv)[ie+1]); 588 } 589 G4double pbval = (xc<=1) ? 590 GetMomentumIntegral(tempData,xc,-1) : 591 GetMomentumIntegral(tempData,1.0,-1); 592 thePBvec->PutValues(ie,theEGrid[ie],pbval); 593 delete[] tempData; 594 } 595 596 fSamplingTable->insert(std::make_pair(theKey,thePhysicsTable)); 597 fPBcut->insert(std::make_pair(theKey,thePBvec)); 598 return; 599 } 600 601 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 602 603 G4double G4PenelopeBremsstrahlungFS::SampleGammaEnergy(G4double energy,const G4Material* mat, 604 const G4double cut) const 605 { 606 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut); 607 if (!(fSamplingTable->count(theKey)) || !(fPBcut->count(theKey))) 608 { 609 G4ExceptionDescription ed; 610 ed << "Unable to retrieve the SamplingTable: " << 611 fSamplingTable->count(theKey) << " " << 612 fPBcut->count(theKey) << G4endl; 613 G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()", 614 "em2014",FatalException,ed); 615 } 616 const G4PhysicsTable* theTableInte = fSamplingTable->find(theKey)->second; 617 const G4PhysicsTable* theTableRed = fReducedXSTable->find(theKey)->second; 618 619 //Find the energy bin using bi-partition 620 std::size_t eBin = 0; 621 G4bool firstOrLastBin = false; 622 623 if (energy < theEGrid[0]) //below first bin 624 { 625 eBin = 0; 626 firstOrLastBin = true; 627 } 628 else if (energy > theEGrid[fNBinsE-1]) //after last bin 629 { 630 eBin = fNBinsE-1; 631 firstOrLastBin = true; 632 } 633 else 634 { 635 std::size_t i=0; 636 std::size_t j=fNBinsE-1; 637 while ((j-i)>1) 638 { 639 std::size_t k = (i+j)/2; 640 if (energy > theEGrid[k]) 641 i = k; 642 else 643 j = k; 644 } 645 eBin = i; 646 } 647 648 //Get the appropriate physics vector 649 const G4PhysicsFreeVector* theVec1 = (G4PhysicsFreeVector*) (*theTableInte)[eBin]; 650 651 //Use a "temporary" vector which contains the linear interpolation of the x spectra 652 //in energy. The temporary vector is thread-local, so that there is no conflict. 653 //This is achieved via G4Cache. The theTempVect is allocated only once per thread 654 //(member variable), but it is overwritten at every call of this method 655 //(because the interpolation factors change!) 656 G4PhysicsFreeVector* theTempVec = fCache.Get(); 657 if (!theTempVec) //First time this thread gets the cache 658 { 659 theTempVec = new G4PhysicsFreeVector(fNBinsX); 660 fCache.Put(theTempVec); 661 // The G4AutoDelete takes care here to clean up the vectors 662 G4AutoDelete::Register(theTempVec); 663 if (fVerbosity > 4) 664 G4cout << "Creating new instance of G4PhysicsFreeVector() on the worker" << G4endl; 665 } 666 667 //theTempVect is allocated only once (member variable), but it is overwritten at 668 //every call of this method (because the interpolation factors change!) 669 if (!firstOrLastBin) 670 { 671 const G4PhysicsFreeVector* theVec2 = (G4PhysicsFreeVector*) (*theTableInte)[eBin+1]; 672 for (std::size_t iloop=0;iloop<fNBinsX;++iloop) 673 { 674 G4double val = (*theVec1)[iloop]+(((*theVec2)[iloop]-(*theVec1)[iloop]))* 675 (energy-theEGrid[eBin])/(theEGrid[eBin+1]-theEGrid[eBin]); 676 theTempVec->PutValues(iloop,theXGrid[iloop],val); 677 } 678 } 679 else //first or last bin, no interpolation 680 { 681 for (std::size_t iloop=0;iloop<fNBinsX;++iloop) 682 theTempVec->PutValues(iloop,theXGrid[iloop],(*theVec1)[iloop]); 683 } 684 685 //Start the game 686 G4double pbcut = (*(fPBcut->find(theKey)->second))[eBin]; 687 688 if (!firstOrLastBin) //linear interpolation on pbcut as well 689 { 690 pbcut = (*(fPBcut->find(theKey)->second))[eBin] + 691 ((*(fPBcut->find(theKey)->second))[eBin+1]-(*(fPBcut->find(theKey)->second))[eBin])* 692 (energy-theEGrid[eBin])/(theEGrid[eBin+1]-theEGrid[eBin]); 693 } 694 695 G4double pCumulative = (*theTempVec)[fNBinsX-1]; //last value 696 697 G4double eGamma = 0; 698 G4int nIterations = 0; 699 do 700 { 701 G4double pt = pbcut + G4UniformRand()*(pCumulative - pbcut); 702 nIterations++; 703 704 //find where it is 705 std::size_t ibin = 0; 706 if (pt < (*theTempVec)[0]) 707 ibin = 0; 708 else if (pt > (*theTempVec)[fNBinsX-1]) 709 { 710 //We observed problems due to numerical rounding here (STT). 711 //delta here is a tiny positive number 712 G4double delta = pt-(*theTempVec)[fNBinsX-1]; 713 if (delta < pt*1e-10) // very small! Numerical rounding only 714 { 715 ibin = fNBinsX-2; 716 G4ExceptionDescription ed; 717 ed << "Found that (pt > (*theTempVec)[fNBinsX-1]) with pt = " << pt << 718 " , (*theTempVec)[fNBinsX-1] = " << (*theTempVec)[fNBinsX-1] << " and delta = " << 719 (pt-(*theTempVec)[fNBinsX-1]) << G4endl; 720 ed << "Possible symptom of problem with numerical precision" << G4endl; 721 G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()", 722 "em2015",JustWarning,ed); 723 } 724 else //real problem 725 { 726 G4ExceptionDescription ed; 727 ed << "Crash at (pt > (*theTempVec)[fNBinsX-1]) with pt = " << pt << 728 " , (*theTempVec)[fNBinsX-1]=" << (*theTempVec)[fNBinsX-1] << " and fNBinsX = " << 729 fNBinsX << G4endl; 730 ed << "Material: " << mat->GetName() << ", energy = " << energy/keV << " keV" << 731 G4endl; 732 G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()", 733 "em2015",FatalException,ed); 734 } 735 } 736 else 737 { 738 std::size_t i=0; 739 std::size_t j=fNBinsX-1; 740 while ((j-i)>1) 741 { 742 std::size_t k = (i+j)/2; 743 if (pt > (*theTempVec)[k]) 744 i = k; 745 else 746 j = k; 747 } 748 ibin = i; 749 } 750 751 G4double w1 = theXGrid[ibin]; 752 G4double w2 = theXGrid[ibin+1]; 753 754 const G4PhysicsFreeVector* v1 = (G4PhysicsFreeVector*) (*theTableRed)[ibin]; 755 const G4PhysicsFreeVector* v2 = (G4PhysicsFreeVector*) (*theTableRed)[ibin+1]; 756 //Remember: the table fReducedXSTable has a fake first point in energy 757 //so, it contains one more bin than fNBinsE. 758 G4double pdf1 = G4Exp((*v1)[eBin+1]); 759 G4double pdf2 = G4Exp((*v2)[eBin+1]); 760 G4double deltaW = w2-w1; 761 G4double dpdfb = pdf2-pdf1; 762 G4double B = dpdfb/deltaW; 763 G4double A = pdf1-B*w1; 764 //I already made an interpolation in energy, so I can use the actual value for the 765 //calculation of the wbcut, instead of the grid values (except for the last bin) 766 G4double wbcut = (cut < energy) ? cut/energy : 1.0; 767 if (firstOrLastBin) //this is an particular case: no interpolation available 768 wbcut = (cut < theEGrid[eBin]) ? cut/theEGrid[eBin] : 1.0; 769 770 if (w1 < wbcut) 771 w1 = wbcut; 772 if (w2 < w1) 773 { 774 //This configuration can happen if initially wbcut > w2 > w1. Due to the previous 775 //statement, (w1 = wbcut), it becomes wbcut = w1 > w2. In this case, it is not a 776 //real problem. It becomes a problem if w2 < w1 before the w1 = wbcut statement. Issue 777 //a warning only in this specific case. 778 if (w2 > wbcut) 779 { 780 G4ExceptionDescription ed; 781 ed << "Warning in G4PenelopeBremsstrahlungFS::SampleX()" << G4endl; 782 ed << "Conflicting end-point values: w1=" << w1 << "; w2 = " << w2 << G4endl; 783 ed << "wbcut = " << wbcut << " energy= " << energy/keV << " keV" << G4endl; 784 ed << "cut = " << cut/keV << " keV" << G4endl; 785 G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()","em2015", 786 JustWarning,ed); 787 } 788 return w1*energy; 789 } 790 791 G4double pmax = std::max(A+B*w1,A+B*w2); 792 G4bool loopAgain = false; 793 do 794 { 795 loopAgain = false; 796 eGamma = w1* std::pow((w2/w1),G4UniformRand()); 797 if (G4UniformRand()*pmax > (A+B*eGamma)) 798 loopAgain = true; 799 }while(loopAgain); 800 eGamma *= energy; 801 if (nIterations > 100) //protection against infinite loops 802 return eGamma; 803 }while(eGamma < cut); //repeat if sampled sub-cut! 804 805 return eGamma; 806 } 807