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 // Thermal Neutron Scattering 27 // Koi, Tatsumi (SLAC/SCCS) 28 // 29 // Class Description: 30 // 31 // Final State Generators for a high precision 32 // libraries) description of themal neutron sc 33 // Based on Thermal neutron scattering files 34 // from the evaluated nuclear data files ENDF/ 35 // To be used in your physics list in case you 36 // In this case you want to register an object 37 // the corresponding process. 38 39 // 070625 Fix memory leaking at destructor by 40 // 081201 Fix memory leaking at destructor by 41 // 100729 Add model name in constructor Proble 42 // P. Arce, June-2014 Conversion neutron_hp to 43 // 44 #include "G4ParticleHPThermalScattering.hh" 45 46 #include "G4ElementTable.hh" 47 #include "G4MaterialTable.hh" 48 #include "G4Neutron.hh" 49 #include "G4ParticleHPElastic.hh" 50 #include "G4ParticleHPManager.hh" 51 #include "G4ParticleHPThermalScatteringData.hh 52 #include "G4ParticleHPThermalScatteringNames.h 53 #include "G4SystemOfUnits.hh" 54 #include "G4Threading.hh" 55 56 G4ParticleHPThermalScattering::G4ParticleHPThe 57 : G4HadronicInteraction("NeutronHPThermalSca 58 { 59 theHPElastic = new G4ParticleHPElastic(); 60 61 SetMinEnergy(0. * eV); 62 SetMaxEnergy(4 * eV); 63 theXSection = new G4ParticleHPThermalScatter 64 65 nMaterial = 0; 66 nElement = 0; 67 } 68 69 G4ParticleHPThermalScattering::~G4ParticleHPTh 70 { 71 delete theHPElastic; 72 } 73 74 void G4ParticleHPThermalScattering::clearCurre 75 { 76 if (incoherentFSs != nullptr) { 77 for (auto it = incoherentFSs->cbegin(); it 78 for (auto itt = it->second->cbegin(); it 79 for (auto ittt = itt->second->cbegin() 80 delete *ittt; 81 } 82 delete itt->second; 83 } 84 delete it->second; 85 } 86 } 87 88 if (coherentFSs != nullptr) { 89 for (auto it = coherentFSs->cbegin(); it ! 90 for (auto itt = it->second->cbegin(); it 91 for (auto ittt = itt->second->cbegin() 92 delete *ittt; 93 } 94 delete itt->second; 95 } 96 delete it->second; 97 } 98 } 99 100 if (inelasticFSs != nullptr) { 101 for (auto it = inelasticFSs->cbegin(); it 102 for (auto itt = it->second->cbegin(); it 103 for (auto ittt = itt->second->cbegin() 104 for (auto it4 = (*ittt)->vE_isoAngle 105 { 106 delete *it4; 107 } 108 delete *ittt; 109 } 110 delete itt->second; 111 } 112 delete it->second; 113 } 114 } 115 116 incoherentFSs = nullptr; 117 coherentFSs = nullptr; 118 inelasticFSs = nullptr; 119 } 120 121 void G4ParticleHPThermalScattering::BuildPhysi 122 { 123 buildPhysicsTable(); 124 theHPElastic->BuildPhysicsTable(particle); 125 } 126 127 std::map<G4double, std::vector<std::pair<G4dou 128 G4ParticleHPThermalScattering::readACoherentFS 129 { 130 auto aCoherentFSDATA = new std::map<G4double 131 132 std::istringstream theChannel(std::ios::in); 133 G4ParticleHPManager::GetInstance()->GetDataS 134 135 std::vector<G4double> vBraggE; 136 137 G4int dummy; 138 while (theChannel >> dummy) // MF // Loop c 139 { 140 theChannel >> dummy; // MT 141 G4double temp; 142 theChannel >> temp; 143 auto anBragE_P = new std::vector<std::pair 144 145 G4int n; 146 theChannel >> n; 147 for (G4int i = 0; i < n; ++i) { 148 G4double Ei; 149 G4double Pi; 150 if (aCoherentFSDATA->empty()) { 151 theChannel >> Ei; 152 vBraggE.push_back(Ei); 153 } 154 else { 155 Ei = vBraggE[i]; 156 } 157 theChannel >> Pi; 158 anBragE_P->push_back(new std::pair<G4dou 159 } 160 aCoherentFSDATA->insert( 161 std::pair<G4double, std::vector<std::pai 162 } 163 return aCoherentFSDATA; 164 } 165 166 std::map<G4double, std::vector<E_P_E_isoAng*>* 167 G4ParticleHPThermalScattering::readAnInelastic 168 { 169 auto anT_E_P_E_isoAng = new std::map<G4doubl 170 171 std::istringstream theChannel(std::ios::in); 172 G4ParticleHPManager::GetInstance()->GetDataS 173 174 G4int dummy; 175 while (theChannel >> dummy) // MF // Loop c 176 { 177 theChannel >> dummy; // MT 178 G4double temp; 179 theChannel >> temp; 180 auto vE_P_E_isoAng = new std::vector<E_P_E 181 G4int n; 182 theChannel >> n; 183 for (G4int i = 0; i < n; ++i) { 184 vE_P_E_isoAng->push_back(readAnE_P_E_iso 185 } 186 anT_E_P_E_isoAng->insert(std::pair<G4doubl 187 } 188 189 return anT_E_P_E_isoAng; 190 } 191 192 E_P_E_isoAng* 193 G4ParticleHPThermalScattering::readAnE_P_E_iso 194 { 195 auto aData = new E_P_E_isoAng; 196 197 G4double dummy; 198 G4double energy; 199 G4int nep, nl; 200 *file >> dummy; 201 *file >> energy; 202 aData->energy = energy * eV; 203 *file >> dummy; 204 *file >> dummy; 205 *file >> nep; 206 *file >> nl; 207 aData->n = nep / nl; 208 for (G4int i = 0; i < aData->n; ++i) { 209 G4double prob; 210 auto anE_isoAng = new E_isoAng; 211 aData->vE_isoAngle.push_back(anE_isoAng); 212 *file >> energy; 213 anE_isoAng->energy = energy * eV; 214 anE_isoAng->n = nl - 2; 215 anE_isoAng->isoAngle.resize(anE_isoAng->n) 216 *file >> prob; 217 aData->prob.push_back(prob); 218 // G4cout << "G4ParticleHPThermalScatterin 219 // << " " << aData->prob[ i ] << G4endl; 220 for (G4int j = 0; j < anE_isoAng->n; ++j) 221 G4double x; 222 *file >> x; 223 anE_isoAng->isoAngle[j] = x; 224 } 225 } 226 227 // Calcuate sum_of_provXdEs 228 G4double total = 0; 229 aData->secondary_energy_cdf.push_back(0.); 230 for (G4int i = 0; i < aData->n - 1; ++i) { 231 G4double E_L = aData->vE_isoAngle[i]->ener 232 G4double E_H = aData->vE_isoAngle[i + 1]-> 233 G4double dE = E_H - E_L; 234 G4double pdf = (aData->prob[i] + aData->pr 235 total += (pdf); 236 aData->secondary_energy_cdf.push_back(tota 237 aData->secondary_energy_pdf.push_back(pdf) 238 aData->secondary_energy_value.push_back(E_ 239 } 240 241 aData->sum_of_probXdEs = total; 242 243 // Normalize CDF 244 aData->secondary_energy_cdf_size = (G4int)aD 245 for (G4int i = 0; i < aData->secondary_energ 246 aData->secondary_energy_cdf[i] /= total; 247 } 248 249 return aData; 250 } 251 252 std::map<G4double, std::vector<E_isoAng*>*>* 253 G4ParticleHPThermalScattering::readAnIncoheren 254 { 255 auto T_E = new std::map<G4double, std::vecto 256 257 // std::ifstream theChannel( name.c_str() ); 258 std::istringstream theChannel(std::ios::in); 259 G4ParticleHPManager::GetInstance()->GetDataS 260 261 G4int dummy; 262 while (theChannel >> dummy) // MF // Loop c 263 { 264 theChannel >> dummy; // MT 265 G4double temp; 266 theChannel >> temp; 267 auto vE_isoAng = new std::vector<E_isoAng* 268 G4int n; 269 theChannel >> n; 270 for (G4int i = 0; i < n; i++) 271 vE_isoAng->push_back(readAnE_isoAng(&the 272 T_E->insert(std::pair<G4double, std::vecto 273 } 274 // theChannel.close(); 275 276 return T_E; 277 } 278 279 E_isoAng* G4ParticleHPThermalScattering::readA 280 { 281 auto aData = new E_isoAng; 282 283 G4double dummy; 284 G4double energy; 285 G4int n; 286 *file >> dummy; 287 *file >> energy; 288 *file >> dummy; 289 *file >> dummy; 290 *file >> n; 291 *file >> dummy; 292 aData->energy = energy * eV; 293 aData->n = n - 2; 294 aData->isoAngle.resize(n); 295 296 *file >> dummy; 297 *file >> dummy; 298 for (G4int i = 0; i < aData->n; i++) 299 *file >> aData->isoAngle[i]; 300 301 return aData; 302 } 303 304 G4HadFinalState* G4ParticleHPThermalScattering 305 306 { 307 // Select Element > Reaction > 308 309 const G4Material* theMaterial = aTrack.GetMa 310 G4double aTemp = theMaterial->GetTemperature 311 auto n = (G4int)theMaterial->GetNumberOfElem 312 313 G4bool findThermalElement = false; 314 G4int ielement; 315 const G4Element* theElement = nullptr; 316 for (G4int i = 0; i < n; ++i) { 317 theElement = theMaterial->GetElement(i); 318 // Select target element 319 if (aNucleus.GetZ_asInt() == (G4int)(theEl 320 // Check Applicability of Thermal Scatte 321 if (getTS_ID(nullptr, theElement) != -1) 322 ielement = getTS_ID(nullptr, theElemen 323 findThermalElement = true; 324 break; 325 } 326 if (getTS_ID(theMaterial, theElement) != 327 ielement = getTS_ID(theMaterial, theEl 328 findThermalElement = true; 329 break; 330 } 331 } 332 } 333 334 if (findThermalElement) { 335 // Select Reaction (Inelastic, coherent, 336 const G4ParticleDefinition* pd = aTrack.Ge 337 auto dp = new G4DynamicParticle(pd, aTrack 338 G4double total = theXSection->GetCrossSect 339 G4double inelastic = theXSection->GetInela 340 341 G4double random = G4UniformRand(); 342 if (random <= inelastic / total) { 343 // Inelastic 344 345 std::vector<G4double> v_temp; 346 v_temp.clear(); 347 for (auto it = inelasticFSs->find(ieleme 348 it != inelasticFSs->find(ielement)- 349 { 350 v_temp.push_back(it->first); 351 } 352 353 std::pair<G4double, G4double> tempLH = f 354 // 355 // For T_L aNEP_EPM_TL and T_H aNEP_EPM 356 // 357 std::vector<E_P_E_isoAng*>* vNEP_EPM_TL 358 std::vector<E_P_E_isoAng*>* vNEP_EPM_TH 359 360 if (tempLH.first != 0.0 && tempLH.second 361 vNEP_EPM_TL = inelasticFSs->find(ielem 362 vNEP_EPM_TH = inelasticFSs->find(ielem 363 } 364 else if (tempLH.first == 0.0) { 365 auto itm = inelasticFSs->find(ielement 366 vNEP_EPM_TL = itm->second; 367 ++itm; 368 vNEP_EPM_TH = itm->second; 369 tempLH.first = tempLH.second; 370 tempLH.second = itm->first; 371 } 372 else if (tempLH.second == 0.0) { 373 auto itm = inelasticFSs->find(ielement 374 --itm; 375 vNEP_EPM_TH = itm->second; 376 --itm; 377 vNEP_EPM_TL = itm->second; 378 tempLH.second = tempLH.first; 379 tempLH.first = itm->first; 380 } 381 382 G4double sE = 0., mu = 1.0; 383 384 // New Geant4 method - Stochastic temper 385 // (continuous temperature interpolation 386 std::pair<G4double, G4double> secondaryP 387 G4double rand_temp = G4UniformRand(); 388 if (rand_temp < (aTemp - tempLH.first) / 389 secondaryParam = sample_inelastic_E_mu 390 else 391 secondaryParam = sample_inelastic_E_mu 392 393 sE = secondaryParam.first; 394 mu = secondaryParam.second; 395 396 // set 397 theParticleChange.SetEnergyChange(sE); 398 G4double phi = CLHEP::twopi * G4UniformR 399 G4double sint = std::sqrt(1 - mu * mu); 400 theParticleChange.SetMomentumChange(sint 401 } 402 else if (random 403 <= (inelastic + theXSection->GetC 404 / total) 405 { 406 // Coherent Elastic 407 408 G4double E = aTrack.GetKineticEnergy(); 409 410 // T_L and T_H 411 std::vector<G4double> v_temp; 412 v_temp.clear(); 413 for (auto it = coherentFSs->find(ielemen 414 it != coherentFSs->find(ielement)-> 415 { 416 v_temp.push_back(it->first); 417 } 418 419 // T_L T_H 420 std::pair<G4double, G4double> tempLH = f 421 // 422 // 423 // For T_L anEPM_TL and T_H anEPM_TH 424 // 425 std::vector<std::pair<G4double, G4double 426 std::vector<std::pair<G4double, G4double 427 428 if (tempLH.first != 0.0 && tempLH.second 429 pvE_p_TL = coherentFSs->find(ielement) 430 pvE_p_TH = coherentFSs->find(ielement) 431 } 432 else if (tempLH.first == 0.0) { 433 pvE_p_TL = coherentFSs->find(ielement) 434 pvE_p_TH = coherentFSs->find(ielement) 435 tempLH.first = tempLH.second; 436 tempLH.second = v_temp[1]; 437 } 438 else if (tempLH.second == 0.0) { 439 pvE_p_TH = coherentFSs->find(ielement) 440 auto itv = v_temp.cend(); 441 --itv; 442 --itv; 443 pvE_p_TL = coherentFSs->find(ielement) 444 tempLH.second = tempLH.first; 445 tempLH.first = *itv; 446 } 447 else { 448 // tempLH.first == 0.0 && tempLH.secon 449 throw G4HadronicException( 450 __FILE__, __LINE__, 451 "A problem is found in Thermal Scatt 452 } 453 454 std::vector<G4double> vE_T; 455 std::vector<G4double> vp_T; 456 457 auto n1 = (G4int)pvE_p_TL->size(); 458 459 // New Geant4 method - Stochastic interp 460 std::vector<std::pair<G4double, G4double 461 G4double rand_temp = G4UniformRand(); 462 if (rand_temp < (aTemp - tempLH.first) / 463 pvE_p_T_sampled = pvE_p_TH; 464 else 465 pvE_p_T_sampled = pvE_p_TL; 466 467 // 171005 fix bug, contribution from H.N 468 for (G4int i = 0; i < n1; ++i) { 469 vE_T.push_back((*pvE_p_T_sampled)[i]-> 470 vp_T.push_back((*pvE_p_T_sampled)[i]-> 471 } 472 473 G4int j = 0; 474 for (G4int i = 1; i < n1; ++i) { 475 if (E / eV < vE_T[i]) { 476 j = i - 1; 477 break; 478 } 479 } 480 481 G4double rand_for_mu = G4UniformRand(); 482 483 G4int k = 0; 484 for (G4int i = 0; i <= j; ++i) { 485 G4double Pi = vp_T[i] / vp_T[j]; 486 if (rand_for_mu < Pi) { 487 k = i; 488 break; 489 } 490 } 491 492 G4double Ei = vE_T[k]; 493 494 G4double mu = 1 - 2 * Ei / (E / eV); 495 496 if (mu < -1.0) mu = -1.0; 497 498 theParticleChange.SetEnergyChange(E); 499 G4double phi = CLHEP::twopi * G4UniformR 500 G4double sint = std::sqrt(1 - mu * mu); 501 theParticleChange.SetMomentumChange(sint 502 } 503 else { 504 // InCoherent Elastic 505 506 // T_L and T_H 507 std::vector<G4double> v_temp; 508 v_temp.clear(); 509 for (auto it = incoherentFSs->find(ielem 510 it != incoherentFSs->find(ielement) 511 { 512 v_temp.push_back(it->first); 513 } 514 515 // T_L T_H 516 std::pair<G4double, G4double> tempLH = f 517 518 // 519 // For T_L anEPM_TL and T_H anEPM_TH 520 // 521 522 E_isoAng anEPM_TL_E; 523 E_isoAng anEPM_TH_E; 524 525 if (tempLH.first != 0.0 && tempLH.second 526 // Interpolate TL and TH 527 anEPM_TL_E = create_E_isoAng_from_ener 528 aTrack.GetKineticEnergy(), 529 incoherentFSs->find(ielement)->secon 530 anEPM_TH_E = create_E_isoAng_from_ener 531 aTrack.GetKineticEnergy(), 532 incoherentFSs->find(ielement)->secon 533 } 534 else if (tempLH.first == 0.0) { 535 // Extrapolate T0 and T1 536 anEPM_TL_E = create_E_isoAng_from_ener 537 aTrack.GetKineticEnergy(), 538 incoherentFSs->find(ielement)->secon 539 anEPM_TH_E = create_E_isoAng_from_ener 540 aTrack.GetKineticEnergy(), 541 incoherentFSs->find(ielement)->secon 542 tempLH.first = tempLH.second; 543 tempLH.second = v_temp[1]; 544 } 545 else if (tempLH.second == 0.0) { 546 // Extrapolate Tmax-1 and Tmax 547 anEPM_TH_E = create_E_isoAng_from_ener 548 aTrack.GetKineticEnergy(), 549 incoherentFSs->find(ielement)->secon 550 auto itv = v_temp.cend(); 551 --itv; 552 --itv; 553 anEPM_TL_E = create_E_isoAng_from_ener 554 aTrack.GetKineticEnergy(), incoheren 555 tempLH.second = tempLH.first; 556 tempLH.first = *itv; 557 } 558 559 // E_isoAng for aTemp and aTrack.GetKine 560 G4double mu = 1.0; 561 562 // New Geant4 method - Stochastic interp 563 E_isoAng anEPM_T_E_sampled; 564 G4double rand_temp = G4UniformRand(); 565 if (rand_temp < (aTemp - tempLH.first) / 566 anEPM_T_E_sampled = std::move(anEPM_TH 567 else 568 anEPM_T_E_sampled = std::move(anEPM_TL 569 570 mu = getMu(&anEPM_T_E_sampled); 571 572 // Set Final State 573 theParticleChange.SetEnergyChange(aTrack 574 G4double phi = CLHEP::twopi * G4UniformR 575 G4double sint = std::sqrt(1 - mu * mu); 576 theParticleChange.SetMomentumChange(sint 577 } 578 delete dp; 579 580 return &theParticleChange; 581 } 582 // Not thermal element 583 // Neutron HP will handle 584 return theHPElastic->ApplyYourself(aTrack, a 585 true); / 586 } 587 588 //******************************************** 589 // Geant4 new algorithm 590 //******************************************** 591 592 //-------------------------------------------- 593 // New method added by L. Thulliez 2021 (CEA-S 594 //-------------------------------------------- 595 std::pair<G4double, G4int> 596 G4ParticleHPThermalScattering::sample_inelasti 597 598 { 599 G4int i = 0; 600 G4double sE_value = 0; 601 602 for (; i < anE_P_E_isoAng->secondary_energy_ 603 if (rndm1 >= anE_P_E_isoAng->secondary_ene 604 && rndm1 < anE_P_E_isoAng->secondary_e 605 { 606 G4double sE_value_i = anE_P_E_isoAng->se 607 G4double sE_pdf_i = anE_P_E_isoAng->seco 608 G4double sE_value_i1 = anE_P_E_isoAng->s 609 G4double sE_pdf_i1 = anE_P_E_isoAng->sec 610 611 G4double lambda = 0; 612 G4double alpha = (sE_pdf_i1 - sE_pdf_i) 613 G4double rndm = rndm1; 614 615 if (std::fabs(alpha) < 1E-8) { 616 lambda = rndm2; 617 } 618 else { 619 G4double beta = 2 * sE_pdf_i / (sE_pdf 620 rndm = rndm2; 621 G4double gamma = -rndm; 622 G4double delta = beta * beta - 4 * alp 623 624 if (delta < 0 && std::fabs(delta) < 1. 625 626 lambda = -beta + std::sqrt(delta); 627 lambda = lambda / (2 * alpha); 628 629 if (lambda > 1) 630 lambda = 1; 631 else if (lambda < 0) 632 lambda = 0; 633 } 634 635 sE_value = sE_value_i + lambda * (sE_val 636 637 break; 638 } 639 } 640 641 return std::pair<G4double, G4int>(sE_value, 642 } 643 644 //-------------------------------------------- 645 // New method added by L. Thulliez 2021 (CEA-S 646 //-------------------------------------------- 647 std::pair<G4double, G4double> 648 G4ParticleHPThermalScattering::sample_inelasti 649 650 { 651 // Sample primary energy bin 652 std::map<G4double, G4int> map_energy; 653 map_energy.clear(); 654 std::vector<G4double> v_energy; 655 v_energy.clear(); 656 G4int i = 0; 657 for (auto itv = vNEP_EPM->cbegin(); itv != v 658 v_energy.push_back((*itv)->energy); 659 map_energy.insert(std::pair<G4double, G4in 660 i++; 661 } 662 663 std::pair<G4double, G4double> energyLH = fin 664 665 std::vector<E_P_E_isoAng*> pE_P_E_isoAng_lim 666 667 if (energyLH.first != 0.0 && energyLH.second 668 pE_P_E_isoAng_limit[0] = (*vNEP_EPM)[map_e 669 pE_P_E_isoAng_limit[1] = (*vNEP_EPM)[map_e 670 } 671 else if (energyLH.first == 0.0) { 672 pE_P_E_isoAng_limit[0] = (*vNEP_EPM)[0]; 673 pE_P_E_isoAng_limit[1] = (*vNEP_EPM)[1]; 674 } 675 if (energyLH.second == 0.0) { 676 pE_P_E_isoAng_limit[1] = (*vNEP_EPM).back( 677 auto itv = vNEP_EPM->cend(); 678 --itv; 679 --itv; 680 pE_P_E_isoAng_limit[0] = *itv; 681 } 682 683 // Compute interpolation factor of the incid 684 G4double factor = (energyLH.second - pE) / ( 685 686 if ((energyLH.second - pE) <= 0. && std::fab 687 if ((energyLH.first - pE) >= 0. && std::fabs 688 689 G4double rndm1 = G4UniformRand(); 690 G4double rndm2 = G4UniformRand(); 691 692 // Sample secondary neutron energy 693 std::pair<G4double, G4int> sE_lower = sample 694 std::pair<G4double, G4int> sE_upper = sample 695 G4double sE = factor * sE_lower.first + (1 - 696 sE = sE * eV; 697 698 // Sample cosine knowing the secondary neutr 699 rndm1 = G4UniformRand(); 700 rndm2 = G4UniformRand(); 701 G4double mu_lower = getMu(rndm1, rndm2, pE_P 702 G4double mu_upper = getMu(rndm1, rndm2, pE_P 703 G4double mu = factor * mu_lower + (1 - facto 704 705 return std::pair<G4double, G4double>(sE, mu) 706 } 707 708 //-------------------------------------------- 709 // New method added by L. Thulliez 2021 (CEA-S 710 //-------------------------------------------- 711 G4double G4ParticleHPThermalScattering::getMu( 712 { 713 G4double result = 0.0; 714 715 auto in = G4int(rndm1 * ((*anEPM).n)); 716 717 if (in != 0) { 718 G4double mu_l = (*anEPM).isoAngle[in - 1]; 719 G4double mu_h = (*anEPM).isoAngle[in]; 720 result = (mu_h - mu_l) * (rndm1 * ((*anEPM 721 } 722 else { 723 G4double x = rndm1 * (*anEPM).n; 724 G4double ratio = 0.5; 725 if (x <= ratio) { 726 G4double mu_l = -1; 727 G4double mu_h = (*anEPM).isoAngle[0]; 728 result = (mu_h - mu_l) * rndm2 + mu_l; 729 } 730 else { 731 G4double mu_l = (*anEPM).isoAngle[(*anEP 732 G4double mu_h = 1; 733 result = (mu_h - mu_l) * rndm2 + mu_l; 734 } 735 } 736 737 return result; 738 } 739 740 //******************************************** 741 // Geant4 previous algorithm 742 //******************************************** 743 744 G4double G4ParticleHPThermalScattering::getMu( 745 { 746 G4double random = G4UniformRand(); 747 G4double result = 0.0; 748 749 auto in = G4int(random * ((*anEPM).n)); 750 751 if (in != 0) { 752 G4double mu_l = (*anEPM).isoAngle[in - 1]; 753 G4double mu_h = (*anEPM).isoAngle[in]; 754 result = (mu_h - mu_l) * (random * ((*anEP 755 } 756 else { 757 G4double x = random * (*anEPM).n; 758 // Bugzilla 1971 759 G4double ratio = 0.5; 760 G4double xx = G4UniformRand(); 761 if (x <= ratio) { 762 G4double mu_l = -1; 763 G4double mu_h = (*anEPM).isoAngle[0]; 764 result = (mu_h - mu_l) * xx + mu_l; 765 } 766 else { 767 G4double mu_l = (*anEPM).isoAngle[(*anEP 768 G4double mu_h = 1; 769 result = (mu_h - mu_l) * xx + mu_l; 770 } 771 } 772 return result; 773 } 774 775 std::pair<G4double, G4double> G4ParticleHPTher 776 777 { 778 G4double LL = 0.0; 779 G4double H = 0.0; 780 781 // v->size() == 1 --> LL=H=v(0) 782 if (aVector->size() == 1) { 783 LL = aVector->front(); 784 H = aVector->front(); 785 } 786 else { 787 // 1) temp < v(0) -> LL=0.0 H=v(0) 788 // 2) v(i-1) < temp <= v(i) -> LL=v(i-1) H 789 // 3) v(imax) < temp -> LL=v(imax) H=0.0 790 for (auto it = aVector->cbegin(); it != aV 791 if (x <= *it) { 792 H = *it; 793 if (it != aVector->cbegin()) { 794 // 2) 795 it--; 796 LL = *it; 797 } 798 else { 799 // 1) 800 LL = 0.0; 801 } 802 break; 803 } 804 } 805 // 3) 806 if (H == 0.0) LL = aVector->back(); 807 } 808 809 return std::pair<G4double, G4double>(LL, H); 810 } 811 812 G4double G4ParticleHPThermalScattering::get_li 813 814 815 { 816 G4double y = 0.0; 817 if (High.first - Low.first != 0) { 818 y = (High.second - Low.second) / (High.fir 819 } 820 else { 821 if (High.second == Low.second) { 822 y = High.second; 823 } 824 else { 825 G4cout << "G4ParticleHPThermalScattering 826 } 827 } 828 829 return y; 830 } 831 832 E_isoAng G4ParticleHPThermalScattering::create 833 834 { 835 E_isoAng anEPM_T_E; 836 837 std::vector<G4double> v_e; 838 v_e.clear(); 839 for (auto iv = vEPM->cbegin(); iv != vEPM->c 840 v_e.push_back((*iv)->energy); 841 842 std::pair<G4double, G4double> energyLH = fin 843 // G4cout << " " << energy/eV << " " << ener 844 845 E_isoAng* panEPM_T_EL = nullptr; 846 E_isoAng* panEPM_T_EH = nullptr; 847 848 if (energyLH.first != 0.0 && energyLH.second 849 for (auto iv = vEPM->cbegin(); iv != vEPM- 850 if (energyLH.first == (*iv)->energy) { 851 panEPM_T_EL = *iv; 852 ++iv; 853 panEPM_T_EH = *iv; 854 break; 855 } 856 } 857 } 858 else if (energyLH.first == 0.0) { 859 panEPM_T_EL = (*vEPM)[0]; 860 panEPM_T_EH = (*vEPM)[1]; 861 } 862 else if (energyLH.second == 0.0) { 863 panEPM_T_EH = (*vEPM).back(); 864 auto iv = vEPM->cend(); 865 --iv; 866 --iv; 867 panEPM_T_EL = *iv; 868 } 869 870 if (panEPM_T_EL != nullptr && panEPM_T_EH != 871 // checking isoAng has proper values or no 872 // Inelastic/FS, the first and last entri 873 if (!(check_E_isoAng(panEPM_T_EL))) panEPM 874 if (!(check_E_isoAng(panEPM_T_EH))) panEPM 875 876 if (panEPM_T_EL->n == panEPM_T_EH->n) { 877 anEPM_T_E.energy = energy; 878 anEPM_T_E.n = panEPM_T_EL->n; 879 880 for (G4int i = 0; i < panEPM_T_EL->n; ++ 881 G4double angle; 882 angle = get_linear_interpolated( 883 energy, std::pair<G4double, G4double 884 std::pair<G4double, G4double>(energy 885 anEPM_T_E.isoAngle.push_back(angle); 886 } 887 } 888 else { 889 G4Exception("G4ParticleHPThermalScatteri 890 JustWarning, 891 "G4ParticleHPThermalScatteri 892 } 893 } 894 else { 895 G4Exception("G4ParticleHPThermalScattering 896 FatalException, "Pointer panEP 897 } 898 899 return anEPM_T_E; 900 } 901 902 G4double 903 G4ParticleHPThermalScattering::get_secondary_e 904 905 { 906 G4double secondary_energy = 0.0; 907 908 G4int n = anE_P_E_isoAng->n; 909 G4double sum_p = 0.0; // sum_p_H 910 G4double sum_p_L = 0.0; 911 912 G4double total = 0.0; 913 914 /* 915 delete for speed up 916 for ( G4int i = 0 ; i < n-1 ; ++i ) 917 { 918 G4double E_L = anE_P_E_isoAng->vE_isoA 919 G4double E_H = anE_P_E_isoAng->vE_isoA 920 G4double dE = E_H - E_L; 921 total += ( ( anE_P_E_isoAng->prob[i] ) 922 } 923 924 if ( std::abs( total - anE_P_E_isoAng->su 925 anE_P_E_isoAng->sum_of_probXdEs << G4endl 926 */ 927 total = anE_P_E_isoAng->sum_of_probXdEs; 928 929 for (G4int i = 0; i < n - 1; ++i) { 930 G4double E_L = anE_P_E_isoAng->vE_isoAngle 931 G4double E_H = anE_P_E_isoAng->vE_isoAngle 932 G4double dE = E_H - E_L; 933 sum_p += ((anE_P_E_isoAng->prob[i]) * dE); 934 935 if (random <= sum_p / total) { 936 secondary_energy = 937 get_linear_interpolated(random, std::p 938 std::pair<G4do 939 secondary_energy = secondary_energy * eV 940 break; 941 } 942 sum_p_L = sum_p; 943 } 944 945 return secondary_energy; 946 } 947 948 std::pair<G4double, E_isoAng> 949 G4ParticleHPThermalScattering::create_sE_and_E 950 G4double rand_for_sE, G4double pE, std::vect 951 { 952 std::map<G4double, G4int> map_energy; 953 map_energy.clear(); 954 std::vector<G4double> v_energy; 955 v_energy.clear(); 956 G4int i = 0; 957 for (auto itv = vNEP_EPM->cbegin(); itv != v 958 v_energy.push_back((*itv)->energy); 959 map_energy.insert(std::pair<G4double, G4in 960 i++; 961 } 962 963 std::pair<G4double, G4double> energyLH = fin 964 965 E_P_E_isoAng* pE_P_E_isoAng_EL = nullptr; 966 E_P_E_isoAng* pE_P_E_isoAng_EH = nullptr; 967 968 if (energyLH.first != 0.0 && energyLH.second 969 pE_P_E_isoAng_EL = (*vNEP_EPM)[map_energy. 970 pE_P_E_isoAng_EH = (*vNEP_EPM)[map_energy. 971 } 972 else if (energyLH.first == 0.0) { 973 pE_P_E_isoAng_EL = (*vNEP_EPM)[0]; 974 pE_P_E_isoAng_EH = (*vNEP_EPM)[1]; 975 } 976 if (energyLH.second == 0.0) { 977 pE_P_E_isoAng_EH = (*vNEP_EPM).back(); 978 auto itv = vNEP_EPM->cend(); 979 --itv; 980 --itv; 981 pE_P_E_isoAng_EL = *itv; 982 } 983 984 G4double sE; 985 G4double sE_L; 986 G4double sE_H; 987 988 sE_L = get_secondary_energy_from_E_P_E_isoAn 989 sE_H = get_secondary_energy_from_E_P_E_isoAn 990 991 sE = get_linear_interpolated(pE, std::pair<G 992 std::pair<G4dou 993 994 E_isoAng E_isoAng_L = create_E_isoAng_from_e 995 E_isoAng E_isoAng_H = create_E_isoAng_from_e 996 997 E_isoAng anE_isoAng; 998 // For defeating warning message from compil 999 anE_isoAng.n = 1; 1000 anE_isoAng.energy = sE; // never used 1001 if (E_isoAng_L.n == E_isoAng_H.n) { 1002 anE_isoAng.n = E_isoAng_L.n; 1003 for (G4int j = 0; j < anE_isoAng.n; ++j) 1004 G4double angle; 1005 angle = 1006 get_linear_interpolated(sE, std::pair 1007 std::pair<G4d 1008 anE_isoAng.isoAngle.push_back(angle); 1009 } 1010 } 1011 else { 1012 throw G4HadronicException(__FILE__, __LIN 1013 } 1014 1015 return std::pair<G4double, E_isoAng>(sE, an 1016 } 1017 1018 void G4ParticleHPThermalScattering::buildPhys 1019 { 1020 // Is rebuild of physics table a necessity 1021 if (nMaterial == G4Material::GetMaterialTab 1022 && nElement == G4Element::GetElementTab 1023 { 1024 return; 1025 } 1026 nMaterial = G4Material::GetMaterialTable()- 1027 nElement = G4Element::GetElementTable()->si 1028 1029 dic.clear(); 1030 std::map<G4String, G4int> co_dic; 1031 1032 // Searching Nist Materials 1033 static G4ThreadLocal G4MaterialTable* theMa 1034 if (theMaterialTable == nullptr) theMateria 1035 std::size_t numberOfMaterials = G4Material: 1036 for (std::size_t i = 0; i < numberOfMateria 1037 G4Material* material = (*theMaterialTable 1038 auto numberOfElements = (G4int)material-> 1039 for (G4int j = 0; j < numberOfElements; + 1040 const G4Element* element = material->Ge 1041 if (names.IsThisThermalElement(material 1042 G4int ts_ID_of_this_geometry; 1043 G4String ts_ndl_name = names.GetTS_ND 1044 if (co_dic.find(ts_ndl_name) != co_di 1045 ts_ID_of_this_geometry = co_dic.fin 1046 } 1047 else { 1048 ts_ID_of_this_geometry = (G4int)co_ 1049 co_dic.insert(std::pair<G4String, G 1050 } 1051 1052 // G4cout << "Neutron HP Thermal Scat 1053 // << material->GetName() << " 1054 // << " as internal thermal sc 1055 // G4endl; 1056 1057 dic.insert(std::pair<std::pair<G4Mate 1058 std::pair<G4Material*, const G4Elem 1059 } 1060 } 1061 } 1062 1063 // Searching TS Elements 1064 auto theElementTable = G4Element::GetElemen 1065 std::size_t numberOfElements = G4Element::G 1066 for (std::size_t i = 0; i < numberOfElement 1067 const G4Element* element = (*theElementTa 1068 if (names.IsThisThermalElement(element->G 1069 if (names.IsThisThermalElement(element- 1070 G4int ts_ID_of_this_geometry; 1071 G4String ts_ndl_name = names.GetTS_ND 1072 if (co_dic.find(ts_ndl_name) != co_di 1073 ts_ID_of_this_geometry = co_dic.fin 1074 } 1075 else { 1076 ts_ID_of_this_geometry = (G4int)co_ 1077 co_dic.insert(std::pair<G4String, G 1078 } 1079 dic.insert(std::pair<std::pair<const 1080 std::pair<const G4Material*, const 1081 ts_ID_of_this_geometry)); 1082 } 1083 } 1084 } 1085 1086 G4cout << G4endl; 1087 G4cout 1088 << "Neutron HP Thermal Scattering: Follow 1089 << G4endl; 1090 for (const auto& it : dic) { 1091 if (it.first.first != nullptr) { 1092 G4cout << "Material " << it.first.first 1093 << it.first.second->GetName() << 1094 << G4endl; 1095 } 1096 else { 1097 G4cout << "Element " << it.first.second 1098 << it.second << G4endl; 1099 } 1100 } 1101 G4cout << G4endl; 1102 1103 // Read Cross Section Data files 1104 1105 G4ParticleHPManager* hpmanager = G4Particle 1106 coherentFSs = hpmanager->GetThermalScatteri 1107 incoherentFSs = hpmanager->GetThermalScatte 1108 inelasticFSs = hpmanager->GetThermalScatter 1109 1110 if (G4Threading::IsMasterThread()) { 1111 clearCurrentFSData(); 1112 1113 if (coherentFSs == nullptr) 1114 coherentFSs = 1115 new std::map<G4int, std::map<G4double 1116 if (incoherentFSs == nullptr) 1117 incoherentFSs = new std::map<G4int, std 1118 if (inelasticFSs == nullptr) 1119 inelasticFSs = new std::map<G4int, std: 1120 1121 G4String dirName; 1122 if (G4FindDataDir("G4NEUTRONHPDATA") == n 1123 throw G4HadronicException( 1124 __FILE__, __LINE__, 1125 "Please setenv G4NEUTRONHPDATA to poi 1126 dirName = G4FindDataDir("G4NEUTRONHPDATA" 1127 1128 for (const auto& it : co_dic) { 1129 G4String tsndlName = it.first; 1130 G4int ts_ID = it.second; 1131 1132 // Coherent 1133 G4String fsName = "/ThermalScattering/C 1134 G4String fileName = dirName + fsName + 1135 coherentFSs->insert( 1136 std::pair<G4int, std::map<G4double, s 1137 ts_ID, readACoherentFSDATA(fileName 1138 1139 // incoherent elastic 1140 fsName = "/ThermalScattering/Incoherent 1141 fileName = dirName + fsName + tsndlName 1142 incoherentFSs->insert(std::pair<G4int, 1143 ts_ID, readAnIncoherentFSDATA(fileNam 1144 1145 // inelastic 1146 fsName = "/ThermalScattering/Inelastic/ 1147 fileName = dirName + fsName + tsndlName 1148 inelasticFSs->insert(std::pair<G4int, s 1149 ts_ID, readAnInelasticFSDATA(fileName 1150 } 1151 1152 hpmanager->RegisterThermalScatteringCoher 1153 hpmanager->RegisterThermalScatteringIncoh 1154 hpmanager->RegisterThermalScatteringInela 1155 } 1156 1157 theXSection->BuildPhysicsTable(*(G4Neutron: 1158 } 1159 1160 G4int G4ParticleHPThermalScattering::getTS_ID 1161 { 1162 G4int result = -1; 1163 if (dic.find(std::pair<const G4Material*, c 1164 result = dic.find(std::pair<const G4Mater 1165 return result; 1166 } 1167 1168 const std::pair<G4double, G4double> G4Particl 1169 { 1170 // max energy non-conservation is mass of h 1171 return std::pair<G4double, G4double>(10.0 * 1172 } 1173 1174 void G4ParticleHPThermalScattering::AddUserTh 1175 1176 { 1177 names.AddThermalElement(nameG4Element, file 1178 theXSection->AddUserThermalScatteringFile(n 1179 buildPhysicsTable(); 1180 } 1181 1182 G4bool G4ParticleHPThermalScattering::check_E 1183 { 1184 G4bool result = false; 1185 1186 G4int n = anE_IsoAng->n; 1187 G4double sum = 0.0; 1188 for (G4int i = 0; i < n; ++i) { 1189 sum += anE_IsoAng->isoAngle[i]; 1190 } 1191 if (sum != 0.0) result = true; 1192 1193 return result; 1194 } 1195 1196 void G4ParticleHPThermalScattering::ModelDesc 1197 { 1198 outFile << "High Precision model based on t 1199 << "evaluated nuclear data librarie 1200 << "on specific materials\n"; 1201 } 1202