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 // G4ParticleHPThermalScatteringData 27 // 28 // 15-Nov-06 First implementation is done by T 29 // 070625 implement clearCurrentXSData to fix 30 // P. Arce, June-2014 Conversion neutron_hp to 31 // ------------------------------------------- 32 33 #include "G4ParticleHPThermalScatteringData.hh 34 35 #include "G4ElementTable.hh" 36 #include "G4Neutron.hh" 37 #include "G4ParticleHPManager.hh" 38 #include "G4SystemOfUnits.hh" 39 #include "G4Threading.hh" 40 41 #include <algorithm> 42 #include <list> 43 44 G4ParticleHPThermalScatteringData::G4ParticleH 45 : G4VCrossSectionDataSet("NeutronHPThermalSc 46 { 47 // Upper limit of neutron energy 48 emax = 4 * eV; 49 SetMinKinEnergy(0 * MeV); 50 SetMaxKinEnergy(emax); 51 52 ke_cache = 0.0; 53 xs_cache = 0.0; 54 element_cache = nullptr; 55 material_cache = nullptr; 56 57 indexOfThermalElement.clear(); 58 59 names = new G4ParticleHPThermalScatteringNam 60 } 61 62 G4ParticleHPThermalScatteringData::~G4Particle 63 { 64 clearCurrentXSData(); 65 66 delete names; 67 } 68 69 G4bool G4ParticleHPThermalScatteringData::IsIs 70 71 72 { 73 G4double eKin = dp->GetKineticEnergy(); 74 if (eKin > 4.0 * eV // GetMaxKinEnergy() 75 || eKin < 0 // GetMinKinEnergy() 76 || dp->GetDefinition() != G4Neutron::Neu 77 return false; 78 79 if (dic.find(std::pair<const G4Material*, co 80 != dic.end() 81 || dic.find(std::pair<const G4Material*, 82 return true; 83 84 return false; 85 } 86 87 G4double G4ParticleHPThermalScatteringData::Ge 88 89 90 91 92 { 93 ke_cache = dp->GetKineticEnergy(); 94 element_cache = element; 95 material_cache = material; 96 G4double xs = GetCrossSection(dp, element, m 97 xs_cache = xs; 98 return xs; 99 } 100 101 void G4ParticleHPThermalScatteringData::clearC 102 { 103 if (coherent != nullptr) { 104 for (auto it = coherent->cbegin(); it != c 105 if (it->second != nullptr) { 106 for (auto itt = it->second->cbegin(); 107 delete itt->second; 108 } 109 } 110 delete it->second; 111 } 112 coherent->clear(); 113 } 114 115 if (incoherent != nullptr) { 116 for (auto it = incoherent->cbegin(); it != 117 if (it->second != nullptr) { 118 for (auto itt = it->second->cbegin(); 119 delete itt->second; 120 } 121 } 122 delete it->second; 123 } 124 incoherent->clear(); 125 } 126 127 if (inelastic != nullptr) { 128 for (auto it = inelastic->cbegin(); it != 129 if (it->second != nullptr) { 130 for (auto itt = it->second->cbegin(); 131 delete itt->second; 132 } 133 } 134 delete it->second; 135 } 136 inelastic->clear(); 137 } 138 } 139 140 G4bool G4ParticleHPThermalScatteringData::IsAp 141 142 { 143 G4bool result = false; 144 145 G4double eKin = aP->GetKineticEnergy(); 146 // Check energy 147 if (eKin < emax) { 148 // Check Particle Species 149 if (aP->GetDefinition() == G4Neutron::Neut 150 // anEle is one of Thermal elements 151 auto ie = (G4int)anEle->GetIndex(); 152 for (int it : indexOfThermalElement) { 153 if (ie == it) return true; 154 } 155 } 156 } 157 158 return result; 159 } 160 161 void G4ParticleHPThermalScatteringData::BuildP 162 { 163 if (&aP != G4Neutron::Neutron()) 164 throw G4HadronicException(__FILE__, __LINE 165 "Attempt to use 166 167 // std::map < std::pair < G4Material* , cons 168 // 169 dic.clear(); 170 if (G4Threading::IsMasterThread()) clearCurr 171 172 std::map<G4String, G4int> co_dic; 173 174 // Searching Nist Materials 175 static G4ThreadLocal G4MaterialTable* theMat 176 if (theMaterialTable == nullptr) theMaterial 177 std::size_t numberOfMaterials = G4Material:: 178 for (std::size_t i = 0; i < numberOfMaterial 179 G4Material* material = (*theMaterialTable) 180 auto numberOfElements = (G4int)material->G 181 for (G4int j = 0; j < numberOfElements; ++ 182 const G4Element* element = material->Get 183 if (names->IsThisThermalElement(material 184 G4int ts_ID_of_this_geometry; 185 G4String ts_ndl_name = names->GetTS_ND 186 if (co_dic.find(ts_ndl_name) != co_dic 187 ts_ID_of_this_geometry = co_dic.find 188 } 189 else { 190 ts_ID_of_this_geometry = (G4int)co_d 191 co_dic.insert(std::pair<G4String, G4 192 } 193 194 dic.insert(std::pair<std::pair<G4Mater 195 std::pair<G4Material*, const G4Eleme 196 } 197 } 198 } 199 200 // Searching TS Elements 201 auto theElementTable = G4Element::GetElement 202 std::size_t numberOfElements = G4Element::Ge 203 204 for (std::size_t i = 0; i < numberOfElements 205 const G4Element* element = (*theElementTab 206 if (names->IsThisThermalElement(element->G 207 if (names->IsThisThermalElement(element- 208 G4int ts_ID_of_this_geometry; 209 G4String ts_ndl_name = names->GetTS_ND 210 if (co_dic.find(ts_ndl_name) != co_dic 211 ts_ID_of_this_geometry = co_dic.find 212 } 213 else { 214 ts_ID_of_this_geometry = (G4int)co_d 215 co_dic.insert(std::pair<G4String, G4 216 } 217 218 dic.insert(std::pair<std::pair<const G 219 std::pair<const G4Material*, const G 220 ts_ID_of_this_geometry)); 221 } 222 } 223 } 224 225 G4cout << G4endl; 226 G4cout << "Neutron HP Thermal Scattering Dat 227 "are registered." 228 << G4endl; 229 for (const auto& it : dic) { 230 if (it.first.first != nullptr) { 231 G4cout << "Material " << it.first.first- 232 << it.first.second->GetName() << 233 << G4endl; 234 } 235 else { 236 G4cout << "Element " << it.first.second- 237 << it.second << G4endl; 238 } 239 } 240 G4cout << G4endl; 241 242 G4ParticleHPManager* hpmanager = G4ParticleH 243 244 coherent = hpmanager->GetThermalScatteringCo 245 incoherent = hpmanager->GetThermalScattering 246 inelastic = hpmanager->GetThermalScatteringI 247 248 if (G4Threading::IsMasterThread()) { 249 if (coherent == nullptr) 250 coherent = new std::map<G4int, std::map< 251 if (incoherent == nullptr) 252 incoherent = new std::map<G4int, std::ma 253 if (inelastic == nullptr) 254 inelastic = new std::map<G4int, std::map 255 256 // Read Cross Section Data files 257 258 G4String dirName; 259 if (G4FindDataDir("G4NEUTRONHPDATA") == nu 260 throw G4HadronicException( 261 __FILE__, __LINE__, 262 "Please setenv G4NEUTRONHPDATA to poin 263 G4String baseName = G4FindDataDir("G4NEUTR 264 265 dirName = baseName + "/ThermalScattering"; 266 267 G4String ndl_filename; 268 G4String full_name; 269 270 for (const auto& it : co_dic) { 271 ndl_filename = it.first; 272 G4int ts_ID = it.second; 273 274 // Coherent 275 full_name = dirName + "/Coherent/CrossSe 276 auto coh_amapTemp_EnergyCross = readData 277 coherent->insert(std::pair<G4int, std::m 278 ts_ID, coh_amapTemp_EnergyCross)); 279 280 // Incoherent 281 full_name = dirName + "/Incoherent/Cross 282 auto incoh_amapTemp_EnergyCross = readDa 283 incoherent->insert(std::pair<G4int, std: 284 ts_ID, incoh_amapTemp_EnergyCross)); 285 286 // Inelastic 287 full_name = dirName + "/Inelastic/CrossS 288 auto inela_amapTemp_EnergyCross = readDa 289 inelastic->insert(std::pair<G4int, std:: 290 ts_ID, inela_amapTemp_EnergyCross)); 291 } 292 hpmanager->RegisterThermalScatteringCohere 293 hpmanager->RegisterThermalScatteringIncohe 294 hpmanager->RegisterThermalScatteringInelas 295 } 296 } 297 298 std::map<G4double, G4ParticleHPVector*>* 299 G4ParticleHPThermalScatteringData::readData(co 300 { 301 auto aData = new std::map<G4double, G4Partic 302 303 std::istringstream theChannel; 304 G4ParticleHPManager::GetInstance()->GetDataS 305 306 G4int dummy; 307 while (theChannel >> dummy) // MF // Loop c 308 { 309 theChannel >> dummy; // MT 310 G4double temp; 311 theChannel >> temp; 312 auto anEnergyCross = new G4ParticleHPVecto 313 G4int nData; 314 theChannel >> nData; 315 anEnergyCross->Init(theChannel, nData, eV, 316 aData->insert(std::pair<G4double, G4Partic 317 } 318 319 return aData; 320 } 321 322 void G4ParticleHPThermalScatteringData::DumpPh 323 { 324 if (&aP != G4Neutron::Neutron()) 325 throw G4HadronicException(__FILE__, __LINE 326 "Attempt to use 327 } 328 329 G4double G4ParticleHPThermalScatteringData::Ge 330 331 332 { 333 G4double result = 0; 334 335 G4int ts_id = getTS_ID(aM, anE); 336 337 if (ts_id == -1) return result; 338 339 G4double aT = aM->GetTemperature(); 340 341 G4double Xcoh = GetX(aP, aT, coherent->find( 342 G4double Xincoh = GetX(aP, aT, incoherent->f 343 G4double Xinela = GetX(aP, aT, inelastic->fi 344 345 result = Xcoh + Xincoh + Xinela; 346 347 return result; 348 } 349 350 G4double G4ParticleHPThermalScatteringData::Ge 351 352 353 { 354 G4double result = 0; 355 G4int ts_id = getTS_ID(aM, anE); 356 G4double aT = aM->GetTemperature(); 357 result = GetX(aP, aT, inelastic->find(ts_id) 358 return result; 359 } 360 361 G4double G4ParticleHPThermalScatteringData::Ge 362 363 364 { 365 G4double result = 0; 366 G4int ts_id = getTS_ID(aM, anE); 367 G4double aT = aM->GetTemperature(); 368 result = GetX(aP, aT, coherent->find(ts_id)- 369 return result; 370 } 371 372 G4double G4ParticleHPThermalScatteringData::Ge 373 374 375 { 376 G4double result = 0; 377 G4int ts_id = getTS_ID(aM, anE); 378 G4double aT = aM->GetTemperature(); 379 result = GetX(aP, aT, incoherent->find(ts_id 380 return result; 381 } 382 383 G4int G4ParticleHPThermalScatteringData::getTS 384 385 { 386 G4int result = -1; 387 if (dic.find(std::pair<const G4Material*, co 388 != dic.end()) 389 return dic.find(std::pair<const G4Material 390 ->second; 391 if (dic.find(std::pair<const G4Material*, co 392 return dic.find(std::pair<const G4Material 393 return result; 394 } 395 396 G4double G4ParticleHPThermalScatteringData:: 397 GetX(const G4DynamicParticle* aP, G4double aT, 398 std::map<G4double, G4ParticleHPVector*>* 399 { 400 G4double result = 0; 401 if (amapTemp_EnergyCross->empty()) return re 402 403 G4double eKinetic = aP->GetKineticEnergy(); 404 405 if (amapTemp_EnergyCross->size() == 1) { 406 if (std::fabs(aT - amapTemp_EnergyCross->c 407 > 0.1) 408 { 409 G4cout 410 << "G4ParticleHPThermalScatteringData: 411 << "K) is different more than 10% from 412 << amapTemp_EnergyCross->begin()->firs 413 } 414 result = amapTemp_EnergyCross->begin()->se 415 return result; 416 } 417 418 auto it = amapTemp_EnergyCross->cbegin(); 419 for (it = amapTemp_EnergyCross->cbegin(); it 420 if (aT < it->first) break; 421 } 422 if (it == amapTemp_EnergyCross->cbegin()) { 423 ++it; // lower than the first 424 } 425 else if (it == amapTemp_EnergyCross->cend()) 426 --it; // upper than the last 427 } 428 429 G4double TH = it->first; 430 G4double XH = it->second->GetXsec(eKinetic); 431 432 if (it != amapTemp_EnergyCross->cbegin()) -- 433 G4double TL = it->first; 434 G4double XL = it->second->GetXsec(eKinetic); 435 436 if (TH == TL) throw G4HadronicException(__FI 437 438 G4double T = aT; 439 G4double X = (XH - XL) / (TH - TL) * (T - TL 440 result = X; 441 442 return result; 443 } 444 445 void G4ParticleHPThermalScatteringData::AddUse 446 447 { 448 names->AddThermalElement(nameG4Element, file 449 } 450 451 void G4ParticleHPThermalScatteringData::CrossS 452 { 453 outFile << "High Precision cross data based 454 "libraries for neutrons below 5eV 455 } 456