| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4ParticleHPThermalScatteringData << 26 // Thermal Neutron Scattering
>> 27 // Koi, Tatsumi (SCCS/SLAC)
27 // 28 //
>> 29 // Class Description
>> 30 // Cross Sections for a high precision (based on evaluated data
>> 31 // libraries) description of themal neutron scattering below 4 eV;
>> 32 // Based on Thermal neutron scattering files
>> 33 // from the evaluated nuclear data files ENDF/B-VI, Release2
>> 34 // To be used in your physics list in case you need this physics.
>> 35 // In this case you want to register an object of this class with
>> 36 // the corresponding process.
>> 37 // Class Description - End
>> 38
28 // 15-Nov-06 First implementation is done by T 39 // 15-Nov-06 First implementation is done by T. Koi (SLAC/SCCS)
29 // 070625 implement clearCurrentXSData to fix 40 // 070625 implement clearCurrentXSData to fix memory leaking by T. Koi
30 // P. Arce, June-2014 Conversion neutron_hp to 41 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
31 // ------------------------------------------- << 42 //
32 43
33 #include "G4ParticleHPThermalScatteringData.hh << 44 #include <list>
>> 45 #include <algorithm>
34 46
35 #include "G4ElementTable.hh" << 47 #include "G4ParticleHPThermalScatteringData.hh"
36 #include "G4Neutron.hh" <<
37 #include "G4ParticleHPManager.hh" 48 #include "G4ParticleHPManager.hh"
>> 49
38 #include "G4SystemOfUnits.hh" 50 #include "G4SystemOfUnits.hh"
39 #include "G4Threading.hh" << 51 #include "G4Neutron.hh"
>> 52 #include "G4ElementTable.hh"
40 53
41 #include <algorithm> << 54 #include "G4Threading.hh"
42 #include <list> <<
43 55
44 G4ParticleHPThermalScatteringData::G4ParticleH 56 G4ParticleHPThermalScatteringData::G4ParticleHPThermalScatteringData()
45 : G4VCrossSectionDataSet("NeutronHPThermalSc << 57 :G4VCrossSectionDataSet("NeutronHPThermalScatteringData")
46 { << 58 ,coherent(NULL)
47 // Upper limit of neutron energy << 59 ,incoherent(NULL)
48 emax = 4 * eV; << 60 ,inelastic(NULL)
49 SetMinKinEnergy(0 * MeV); << 61 {
50 SetMaxKinEnergy(emax); << 62 // Upper limit of neutron energy
>> 63 emax = 4*eV;
>> 64 SetMinKinEnergy( 0*MeV );
>> 65 SetMaxKinEnergy( emax );
>> 66
>> 67 ke_cache = 0.0;
>> 68 xs_cache = 0.0;
>> 69 element_cache = NULL;
>> 70 material_cache = NULL;
51 71
52 ke_cache = 0.0; << 72 indexOfThermalElement.clear();
53 xs_cache = 0.0; <<
54 element_cache = nullptr; <<
55 material_cache = nullptr; <<
56 73
57 indexOfThermalElement.clear(); << 74 names = new G4ParticleHPThermalScatteringNames();
58 <<
59 names = new G4ParticleHPThermalScatteringNam <<
60 } 75 }
61 76
62 G4ParticleHPThermalScatteringData::~G4Particle 77 G4ParticleHPThermalScatteringData::~G4ParticleHPThermalScatteringData()
63 { 78 {
64 clearCurrentXSData(); <<
65 79
66 delete names; << 80 clearCurrentXSData();
>> 81
>> 82 delete names;
>> 83 }
>> 84
>> 85 G4bool G4ParticleHPThermalScatteringData::IsIsoApplicable( const G4DynamicParticle* dp ,
>> 86 G4int /*Z*/ , G4int /*A*/ ,
>> 87 const G4Element* element ,
>> 88 const G4Material* material )
>> 89 {
>> 90 G4double eKin = dp->GetKineticEnergy();
>> 91 if ( eKin > 4.0*eV //GetMaxKinEnergy()
>> 92 || eKin < 0 //GetMinKinEnergy()
>> 93 || dp->GetDefinition() != G4Neutron::Neutron() ) return false;
>> 94
>> 95 if ( dic.find( std::pair < const G4Material* , const G4Element* > ( (G4Material*)NULL , element ) ) != dic.end()
>> 96 || dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) ) != dic.end() ) return true;
>> 97
>> 98 return false;
>> 99
>> 100 // return IsApplicable( dp , element );
>> 101 /*
>> 102 G4double eKin = dp->GetKineticEnergy();
>> 103 if ( eKin > 4.0*eV //GetMaxKinEnergy()
>> 104 || eKin < 0 //GetMinKinEnergy()
>> 105 || dp->GetDefinition() != G4Neutron::Neutron() ) return false;
>> 106 return true;
>> 107 */
>> 108 }
>> 109
>> 110 G4double G4ParticleHPThermalScatteringData::GetIsoCrossSection( const G4DynamicParticle* dp ,
>> 111 G4int /*Z*/ , G4int /*A*/ ,
>> 112 const G4Isotope* /*iso*/ ,
>> 113 const G4Element* element ,
>> 114 const G4Material* material )
>> 115 {
>> 116 //if ( dp->GetKineticEnergy() == ke_cache && element == element_cache && material == material_cache ) return xs_cache;
>> 117
>> 118 ke_cache = dp->GetKineticEnergy();
>> 119 element_cache = element;
>> 120 material_cache = material;
>> 121 //G4double xs = GetCrossSection( dp , element , material->GetTemperature() );
>> 122 G4double xs = GetCrossSection( dp , element , material );
>> 123 xs_cache = xs;
>> 124 return xs;
>> 125 //return GetCrossSection( dp , element , material->GetTemperature() );
67 } 126 }
68 127
69 G4bool G4ParticleHPThermalScatteringData::IsIs << 128 void G4ParticleHPThermalScatteringData::clearCurrentXSData()
70 <<
71 <<
72 { 129 {
73 G4double eKin = dp->GetKineticEnergy(); << 130 std::map< G4int , std::map< G4double , G4ParticleHPVector* >* >::iterator it;
74 if (eKin > 4.0 * eV // GetMaxKinEnergy() << 131 std::map< G4double , G4ParticleHPVector* >::iterator itt;
75 || eKin < 0 // GetMinKinEnergy() << 132
76 || dp->GetDefinition() != G4Neutron::Neu << 133 if ( coherent != NULL ) {
77 return false; << 134 for ( it = coherent->begin() ; it != coherent->end() ; it++ )
>> 135 {
>> 136 if ( it->second != NULL )
>> 137 {
>> 138 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
>> 139 {
>> 140 delete itt->second;
>> 141 }
>> 142 }
>> 143 delete it->second;
>> 144 }
>> 145 coherent->clear();
>> 146 }
>> 147
>> 148 if ( incoherent != NULL ) {
>> 149 for ( it = incoherent->begin() ; it != incoherent->end() ; it++ )
>> 150 {
>> 151 if ( it->second != NULL )
>> 152 {
>> 153 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
>> 154 {
>> 155 delete itt->second;
>> 156 }
>> 157 }
>> 158 delete it->second;
>> 159 }
>> 160 incoherent->clear();
>> 161 }
78 162
79 if (dic.find(std::pair<const G4Material*, co << 163 if ( inelastic != NULL ) {
80 != dic.end() << 164 for ( it = inelastic->begin() ; it != inelastic->end() ; it++ )
81 || dic.find(std::pair<const G4Material*, << 165 {
82 return true; << 166 if ( it->second != NULL )
>> 167 {
>> 168 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
>> 169 {
>> 170 delete itt->second;
>> 171 }
>> 172 }
>> 173 delete it->second;
>> 174 }
>> 175 inelastic->clear();
>> 176 }
83 177
84 return false; <<
85 } 178 }
86 179
87 G4double G4ParticleHPThermalScatteringData::Ge << 180
88 << 181
89 << 182 G4bool G4ParticleHPThermalScatteringData::IsApplicable( const G4DynamicParticle* aP , const G4Element* anEle )
90 <<
91 <<
92 { 183 {
93 ke_cache = dp->GetKineticEnergy(); << 184 G4bool result = false;
94 element_cache = element; << 185
95 material_cache = material; << 186 G4double eKin = aP->GetKineticEnergy();
96 G4double xs = GetCrossSection(dp, element, m << 187 // Check energy
97 xs_cache = xs; << 188 if ( eKin < emax )
98 return xs; << 189 {
>> 190 // Check Particle Species
>> 191 if ( aP->GetDefinition() == G4Neutron::Neutron() )
>> 192 {
>> 193 // anEle is one of Thermal elements
>> 194 G4int ie = (G4int) anEle->GetIndex();
>> 195 std::vector < G4int >::iterator it;
>> 196 for ( it = indexOfThermalElement.begin() ; it != indexOfThermalElement.end() ; it++ )
>> 197 {
>> 198 if ( ie == *it ) return true;
>> 199 }
>> 200 }
>> 201 }
>> 202
>> 203 /*
>> 204 if ( names->IsThisThermalElement ( anEle->GetName() ) )
>> 205 {
>> 206 // Check energy and projectile species
>> 207 G4double eKin = aP->GetKineticEnergy();
>> 208 if ( eKin < emax && aP->GetDefinition() == G4Neutron::Neutron() ) result = true;
>> 209 }
>> 210 */
>> 211 return result;
99 } 212 }
100 213
101 void G4ParticleHPThermalScatteringData::clearC << 214
>> 215 void G4ParticleHPThermalScatteringData::BuildPhysicsTable(const G4ParticleDefinition& aP)
102 { 216 {
103 if (coherent != nullptr) { << 217
104 for (auto it = coherent->cbegin(); it != c << 218 if ( &aP != G4Neutron::Neutron() )
105 if (it->second != nullptr) { << 219 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
106 for (auto itt = it->second->cbegin(); << 220
107 delete itt->second; << 221 //std::map < std::pair < G4Material* , const G4Element* > , G4int > dic;
108 } << 222 //
>> 223 dic.clear();
>> 224 if ( G4Threading::IsMasterThread() ) clearCurrentXSData();
>> 225
>> 226 std::map < G4String , G4int > co_dic;
>> 227
>> 228 //Searching Nist Materials
>> 229 static G4ThreadLocal G4MaterialTable* theMaterialTable = 0 ; if (!theMaterialTable) theMaterialTable= G4Material::GetMaterialTable();
>> 230 size_t numberOfMaterials = G4Material::GetNumberOfMaterials();
>> 231 for ( size_t i = 0 ; i < numberOfMaterials ; i++ )
>> 232 {
>> 233 G4Material* material = (*theMaterialTable)[i];
>> 234 size_t numberOfElements = material->GetNumberOfElements();
>> 235 for ( size_t j = 0 ; j < numberOfElements ; j++ )
>> 236 {
>> 237 const G4Element* element = material->GetElement(j);
>> 238 if ( names->IsThisThermalElement ( material->GetName() , element->GetName() ) )
>> 239 {
>> 240 G4int ts_ID_of_this_geometry;
>> 241 G4String ts_ndl_name = names->GetTS_NDL_Name( material->GetName() , element->GetName() );
>> 242 if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
>> 243 {
>> 244 ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) -> second;
>> 245 }
>> 246 else
>> 247 {
>> 248 ts_ID_of_this_geometry = co_dic.size();
>> 249 co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
>> 250 }
>> 251
>> 252 //G4cout << "Neutron HP Thermal Scattering Data : Registering a material-element pair of "
>> 253 // << material->GetName() << " " << element->GetName()
>> 254 // << " as internal thermal scattering id of " << ts_ID_of_this_geometry << "." << G4endl;
>> 255
>> 256 dic.insert( std::pair < std::pair < G4Material* , const G4Element* > , G4int > ( std::pair < G4Material* , const G4Element* > ( material , element ) , ts_ID_of_this_geometry ) );
>> 257 }
109 } 258 }
110 delete it->second; << 259 }
111 } << 260
112 coherent->clear(); << 261 //Searching TS Elements
113 } << 262 static G4ThreadLocal G4ElementTable* theElementTable = 0 ; if (!theElementTable) theElementTable= G4Element::GetElementTable();
114 << 263 size_t numberOfElements = G4Element::GetNumberOfElements();
115 if (incoherent != nullptr) { << 264 //size_t numberOfThermalElements = 0;
116 for (auto it = incoherent->cbegin(); it != << 265 for ( size_t i = 0 ; i < numberOfElements ; i++ )
117 if (it->second != nullptr) { << 266 {
118 for (auto itt = it->second->cbegin(); << 267 const G4Element* element = (*theElementTable)[i];
119 delete itt->second; << 268 if ( names->IsThisThermalElement ( element->GetName() ) )
120 } << 269 {
>> 270 if ( names->IsThisThermalElement ( element->GetName() ) )
>> 271 {
>> 272 G4int ts_ID_of_this_geometry;
>> 273 G4String ts_ndl_name = names->GetTS_NDL_Name( element->GetName() );
>> 274 if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
>> 275 {
>> 276 ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) -> second;
>> 277 }
>> 278 else
>> 279 {
>> 280 ts_ID_of_this_geometry = co_dic.size();
>> 281 co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
>> 282 }
>> 283
>> 284 //G4cout << "Neutron HP Thermal Scattering: Registering an element of "
>> 285 // << material->GetName() << " " << element->GetName()
>> 286 // << " as internal thermal scattering id of " << ts_ID_of_this_geometry << "." << G4endl;
>> 287
>> 288 dic.insert( std::pair < std::pair < const G4Material* , const G4Element* > , G4int > ( std::pair < const G4Material* , const G4Element* > ( (G4Material*)NULL , element ) , ts_ID_of_this_geometry ) );
>> 289 }
121 } 290 }
122 delete it->second; << 291 }
123 } << 292
124 incoherent->clear(); << 293 G4cout << G4endl;
125 } << 294 G4cout << "Neutron HP Thermal Scattering Data: Following material-element pairs and/or elements are registered." << G4endl;
126 << 295 for ( std::map < std::pair < const G4Material* , const G4Element* > , G4int >::iterator it = dic.begin() ; it != dic.end() ; it++ )
127 if (inelastic != nullptr) { << 296 {
128 for (auto it = inelastic->cbegin(); it != << 297 if ( it->first.first != NULL )
129 if (it->second != nullptr) { << 298 {
130 for (auto itt = it->second->cbegin(); << 299 G4cout << "Material " << it->first.first->GetName() << " - Element " << it->first.second->GetName() << ", internal thermal scattering id " << it->second << G4endl;
131 delete itt->second; <<
132 } <<
133 } 300 }
134 delete it->second; << 301 else
135 } << 302 {
136 inelastic->clear(); << 303 G4cout << "Element " << it->first.second->GetName() << ", internal thermal scattering id " << it->second << G4endl;
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 } 304 }
155 } << 305 }
156 } << 306 G4cout << G4endl;
>> 307
>> 308
>> 309 //G4cout << "Neutron HP Thermal Scattering Data: Following NDL thermal scattering files are assigned to the internal thermal scattering id." << G4endl;
>> 310 //for ( std::map < G4String , G4int >::iterator it = co_dic.begin() ; it != co_dic.end() ; it++ )
>> 311 //{
>> 312 // G4cout << "NDL file name " << it->first << ", internal thermal scattering id " << it->second << G4endl;
>> 313 //}
>> 314
>> 315 G4ParticleHPManager* hpmanager = G4ParticleHPManager::GetInstance();
>> 316
>> 317 coherent = hpmanager->GetThermalScatteringCoherentCrossSections();
>> 318 incoherent = hpmanager->GetThermalScatteringIncoherentCrossSections();
>> 319 inelastic = hpmanager->GetThermalScatteringInelasticCrossSections();
>> 320
>> 321 if ( G4Threading::IsMasterThread() ) {
>> 322
>> 323 if ( coherent == NULL ) coherent = new std::map< G4int , std::map< G4double , G4ParticleHPVector* >* >;
>> 324 if ( incoherent == NULL ) incoherent = new std::map< G4int , std::map< G4double , G4ParticleHPVector* >* >;
>> 325 if ( inelastic == NULL ) inelastic = new std::map< G4int , std::map< G4double , G4ParticleHPVector* >* >;
>> 326
>> 327
>> 328 // Read Cross Section Data files
157 329
158 return result; << 330 G4String dirName;
>> 331 if ( !std::getenv( "G4NEUTRONHPDATA" ) )
>> 332 throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
>> 333 G4String baseName = std::getenv( "G4NEUTRONHPDATA" );
>> 334
>> 335 dirName = baseName + "/ThermalScattering";
>> 336
>> 337 G4String ndl_filename;
>> 338 G4String full_name;
>> 339
>> 340 for ( std::map < G4String , G4int >::iterator it = co_dic.begin() ; it != co_dic.end() ; it++ )
>> 341 {
>> 342
>> 343 ndl_filename = it->first;
>> 344 G4int ts_ID = it->second;
>> 345
>> 346 // Coherent
>> 347 full_name = dirName + "/Coherent/CrossSection/" + ndl_filename;
>> 348 std::map< G4double , G4ParticleHPVector* >* coh_amapTemp_EnergyCross = readData( full_name );
>> 349 coherent->insert ( std::pair < G4int , std::map< G4double , G4ParticleHPVector* >* > ( ts_ID , coh_amapTemp_EnergyCross ) );
>> 350
>> 351 // Incoherent
>> 352 full_name = dirName + "/Incoherent/CrossSection/" + ndl_filename;
>> 353 std::map< G4double , G4ParticleHPVector* >* incoh_amapTemp_EnergyCross = readData( full_name );
>> 354 incoherent->insert ( std::pair < G4int , std::map< G4double , G4ParticleHPVector* >* > ( ts_ID , incoh_amapTemp_EnergyCross ) );
>> 355
>> 356 // Inelastic
>> 357 full_name = dirName + "/Inelastic/CrossSection/" + ndl_filename;
>> 358 std::map< G4double , G4ParticleHPVector* >* inela_amapTemp_EnergyCross = readData( full_name );
>> 359 inelastic->insert ( std::pair < G4int , std::map< G4double , G4ParticleHPVector* >* > ( ts_ID , inela_amapTemp_EnergyCross ) );
>> 360
>> 361 }
>> 362 hpmanager->RegisterThermalScatteringCoherentCrossSections( coherent );
>> 363 hpmanager->RegisterThermalScatteringIncoherentCrossSections( incoherent );
>> 364 hpmanager->RegisterThermalScatteringInelasticCrossSections( inelastic );
>> 365 }
159 } 366 }
160 367
161 void G4ParticleHPThermalScatteringData::BuildP << 368
>> 369
>> 370 std::map< G4double , G4ParticleHPVector* >* G4ParticleHPThermalScatteringData::readData ( G4String full_name )
162 { 371 {
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 372
194 dic.insert(std::pair<std::pair<G4Mater << 373 std::map< G4double , G4ParticleHPVector* >* aData = new std::map< G4double , G4ParticleHPVector* >;
195 std::pair<G4Material*, const G4Eleme << 374
196 } << 375 //std::ifstream theChannel( full_name.c_str() );
197 } << 376 std::istringstream theChannel;
198 } << 377 G4ParticleHPManager::GetInstance()->GetDataStream(full_name,theChannel);
199 378
200 // Searching TS Elements << 379 //G4cout << "G4ParticleHPThermalScatteringData " << name << G4endl;
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 380
225 G4cout << G4endl; << 381 G4int dummy;
226 G4cout << "Neutron HP Thermal Scattering Dat << 382 while ( theChannel >> dummy ) // MF // Loop checking, 11.05.2015, T. Koi
227 "are registered." << 383 {
228 << G4endl; << 384 theChannel >> dummy; // MT
229 for (const auto& it : dic) { << 385 G4double temp;
230 if (it.first.first != nullptr) { << 386 theChannel >> temp;
231 G4cout << "Material " << it.first.first- << 387 G4ParticleHPVector* anEnergyCross = new G4ParticleHPVector;
232 << it.first.second->GetName() << << 388 G4int nData;
233 << G4endl; << 389 theChannel >> nData;
234 } << 390 anEnergyCross->Init ( theChannel , nData , eV , barn );
235 else { << 391 aData->insert ( std::pair < G4double , G4ParticleHPVector* > ( temp , anEnergyCross ) );
236 G4cout << "Element " << it.first.second- << 392 }
237 << it.second << G4endl; << 393 //theChannel.close();
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 394
442 return result; << 395 return aData;
>> 396
>> 397 }
>> 398
>> 399
>> 400
>> 401 void G4ParticleHPThermalScatteringData::DumpPhysicsTable( const G4ParticleDefinition& aP )
>> 402 {
>> 403 if( &aP != G4Neutron::Neutron() )
>> 404 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
>> 405 // G4cout << "G4ParticleHPThermalScatteringData::DumpPhysicsTable still to be implemented"<<G4endl;
>> 406 }
>> 407
>> 408 //#include "G4Nucleus.hh"
>> 409 //#include "G4NucleiPropertiesTable.hh"
>> 410 //#include "G4Neutron.hh"
>> 411 //#include "G4Electron.hh"
>> 412
>> 413
>> 414
>> 415 /*
>> 416 G4double G4ParticleHPThermalScatteringData::GetCrossSection( const G4DynamicParticle* aP , const G4Element*anE , G4double aT )
>> 417 {
>> 418
>> 419 G4double result = 0;
>> 420 const G4Material* aM = NULL;
>> 421
>> 422 G4int iele = anE->GetIndex();
>> 423
>> 424 if ( dic.find( std::pair < const G4Material* , const G4Element* > ( (G4Material*)NULL , anE ) ) != dic.end() )
>> 425 {
>> 426 iele = dic.find( std::pair < const G4Material* , const G4Element* > ( (G4Material*)NULL , anE ) )->second;
>> 427 }
>> 428 else if ( dic.find( std::pair < const G4Material* , const G4Element* > ( aM , anE ) ) != dic.end() )
>> 429 {
>> 430 iele = dic.find( std::pair < const G4Material* , const G4Element* > ( aM , anE ) )->second;
>> 431 }
>> 432 else
>> 433 {
>> 434 return result;
>> 435 }
>> 436
>> 437 G4double Xcoh = GetX ( aP , aT , coherent.find(iele)->second );
>> 438 G4double Xincoh = GetX ( aP , aT , incoherent.find(iele)->second );
>> 439 G4double Xinela = GetX ( aP , aT , inelastic.find(iele)->second );
>> 440
>> 441 result = Xcoh + Xincoh + Xinela;
>> 442
>> 443 //G4cout << "G4ParticleHPThermalScatteringData::GetCrossSection Tot= " << result/barn << " Coherent= " << Xcoh/barn << " Incoherent= " << Xincoh/barn << " Inelastic= " << Xinela/barn << G4endl;
>> 444
>> 445 return result;
>> 446
>> 447 }
>> 448 */
>> 449
>> 450 G4double G4ParticleHPThermalScatteringData::GetCrossSection( const G4DynamicParticle* aP , const G4Element*anE , const G4Material* aM )
>> 451 {
>> 452 G4double result = 0;
>> 453
>> 454 G4int ts_id =getTS_ID( aM , anE );
>> 455
>> 456 if ( ts_id == -1 ) return result;
>> 457
>> 458 G4double aT = aM->GetTemperature();
>> 459
>> 460 G4double Xcoh = GetX ( aP , aT , coherent->find(ts_id)->second );
>> 461 G4double Xincoh = GetX ( aP , aT , incoherent->find(ts_id)->second );
>> 462 G4double Xinela = GetX ( aP , aT , inelastic->find(ts_id)->second );
>> 463
>> 464 result = Xcoh + Xincoh + Xinela;
>> 465
>> 466 //G4cout << "G4ParticleHPThermalScatteringData::GetCrossSection Tot= " << result/barn << " Coherent= " << Xcoh/barn << " Incoherent= " << Xincoh/barn << " Inelastic= " << Xinela/barn << G4endl;
>> 467
>> 468 return result;
>> 469 }
>> 470
>> 471
>> 472 G4double G4ParticleHPThermalScatteringData::GetInelasticCrossSection( const G4DynamicParticle* aP , const G4Element*anE , const G4Material* aM )
>> 473 {
>> 474 G4double result = 0;
>> 475 G4int ts_id = getTS_ID( aM , anE );
>> 476 G4double aT = aM->GetTemperature();
>> 477 result = GetX ( aP , aT , inelastic->find( ts_id )->second );
>> 478 return result;
>> 479 }
>> 480
>> 481 G4double G4ParticleHPThermalScatteringData::GetCoherentCrossSection( const G4DynamicParticle* aP , const G4Element*anE , const G4Material* aM )
>> 482 {
>> 483 G4double result = 0;
>> 484 G4int ts_id = getTS_ID( aM , anE );
>> 485 G4double aT = aM->GetTemperature();
>> 486 result = GetX ( aP , aT , coherent->find( ts_id )->second );
>> 487 return result;
443 } 488 }
444 489
445 void G4ParticleHPThermalScatteringData::AddUse << 490 G4double G4ParticleHPThermalScatteringData::GetIncoherentCrossSection( const G4DynamicParticle* aP , const G4Element*anE , const G4Material* aM )
446 << 491 {
>> 492 G4double result = 0;
>> 493 G4int ts_id = getTS_ID( aM , anE );
>> 494 G4double aT = aM->GetTemperature();
>> 495 result = GetX ( aP , aT , incoherent->find( ts_id )->second );
>> 496 return result;
>> 497 }
>> 498
>> 499
>> 500
>> 501 G4int G4ParticleHPThermalScatteringData::getTS_ID ( const G4Material* material , const G4Element* element )
447 { 502 {
448 names->AddThermalElement(nameG4Element, file << 503 G4int result = -1;
>> 504 if ( dic.find( std::pair < const G4Material* , const G4Element* > ( (G4Material*)NULL , element ) ) != dic.end() )
>> 505 return dic.find( std::pair < const G4Material* , const G4Element* > ( (G4Material*)NULL , element ) )->second;
>> 506 if ( dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) ) != dic.end() )
>> 507 return dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) )->second;
>> 508 return result;
449 } 509 }
450 510
>> 511
>> 512
>> 513
>> 514 G4double G4ParticleHPThermalScatteringData::GetX ( const G4DynamicParticle* aP, G4double aT , std::map < G4double , G4ParticleHPVector* >* amapTemp_EnergyCross )
>> 515 {
>> 516
>> 517 G4double result = 0;
>> 518 if ( amapTemp_EnergyCross->size() == 0 ) return result;
>> 519
>> 520
>> 521 G4double eKinetic = aP->GetKineticEnergy();
>> 522
>> 523 if ( amapTemp_EnergyCross->size() == 1 ) {
>> 524 if ( std::fabs ( aT - amapTemp_EnergyCross->begin()->first ) / amapTemp_EnergyCross->begin()->first > 0.1 ) {
>> 525 G4cout << "G4ParticleHPThermalScatteringData:: The temperature of material ("
>> 526 << aT/kelvin << "K) is different more than 10% from temperature of thermal scattering file expected ("
>> 527 << amapTemp_EnergyCross->begin()->first << "K). Result may not be reliable."
>> 528 << G4endl;
>> 529 }
>> 530 result = amapTemp_EnergyCross->begin()->second->GetXsec ( eKinetic );
>> 531 return result;
>> 532 }
>> 533
>> 534 std::map< G4double , G4ParticleHPVector* >::iterator it;
>> 535 for ( it = amapTemp_EnergyCross->begin() ; it != amapTemp_EnergyCross->end() ; it++ ) {
>> 536 if ( aT < it->first ) break;
>> 537 }
>> 538 //if ( it == amapTemp_EnergyCross->begin() && it != amapTemp_EnergyCross->end() ) it++; // lower than the first
>> 539 //if ( it != amapTemp_EnergyCross->begin() && it == amapTemp_EnergyCross->end() ) it--; // upper than the last
>> 540 if ( it == amapTemp_EnergyCross->begin() ) {
>> 541 it++; // lower than the first
>> 542 } else if ( it == amapTemp_EnergyCross->end() ) {
>> 543 it--; // upper than the last
>> 544 }
>> 545
>> 546 G4double TH = it->first;
>> 547 G4double XH = it->second->GetXsec ( eKinetic );
>> 548
>> 549 //G4cout << "G4ParticleHPThermalScatteringData::GetX TH " << TH << " E " << eKinetic << " XH " << XH << G4endl;
>> 550
>> 551 if ( it != amapTemp_EnergyCross->begin() ) it--;
>> 552 G4double TL = it->first;
>> 553 G4double XL = it->second->GetXsec ( eKinetic );
>> 554
>> 555 //G4cout << "G4ParticleHPThermalScatteringData::GetX TL " << TL << " E " << eKinetic << " XL " << XL << G4endl;
>> 556
>> 557 if ( TH == TL )
>> 558 throw G4HadronicException(__FILE__, __LINE__, "Thermal Scattering Data Error!");
>> 559
>> 560 G4double T = aT;
>> 561 G4double X = ( XH - XL ) / ( TH - TL ) * ( T - TL ) + XL;
>> 562 result = X;
>> 563
>> 564 return result;
>> 565 }
>> 566
>> 567
>> 568 void G4ParticleHPThermalScatteringData::AddUserThermalScatteringFile( G4String nameG4Element , G4String filename )
>> 569 {
>> 570 names->AddThermalElement( nameG4Element , filename );
>> 571 }
451 void G4ParticleHPThermalScatteringData::CrossS 572 void G4ParticleHPThermalScatteringData::CrossSectionDescription(std::ostream& outFile) const
452 { 573 {
453 outFile << "High Precision cross data based << 574 outFile << "High Precision cross data based on thermal scattering data in evaluated nuclear data libraries for neutrons below 5eV on specific materials\n" ;
454 "libraries for neutrons below 5eV <<
455 } 575 }
456 576