00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include <list>
00043 #include <algorithm>
00044
00045 #include "G4NeutronHPThermalScatteringData.hh"
00046 #include "G4SystemOfUnits.hh"
00047 #include "G4Neutron.hh"
00048 #include "G4ElementTable.hh"
00049
00050
00051 G4NeutronHPThermalScatteringData::G4NeutronHPThermalScatteringData()
00052 :G4VCrossSectionDataSet("NeutronHPThermalScatteringData")
00053 {
00054
00055 emax = 4*eV;
00056 SetMinKinEnergy( 0*MeV );
00057 SetMaxKinEnergy( emax );
00058
00059 ke_cache = 0.0;
00060 xs_cache = 0.0;
00061 element_cache = NULL;
00062 material_cache = NULL;
00063
00064 indexOfThermalElement.clear();
00065
00066 names = new G4NeutronHPThermalScatteringNames();
00067
00068
00069 }
00070
00071 G4NeutronHPThermalScatteringData::~G4NeutronHPThermalScatteringData()
00072 {
00073
00074 clearCurrentXSData();
00075
00076 delete names;
00077 }
00078
00079 G4bool G4NeutronHPThermalScatteringData::IsIsoApplicable( const G4DynamicParticle* dp ,
00080 G4int , G4int ,
00081 const G4Element* element ,
00082 const G4Material* material )
00083 {
00084 G4double eKin = dp->GetKineticEnergy();
00085 if ( eKin > 4.0*eV
00086 || eKin < 0
00087 || dp->GetDefinition() != G4Neutron::Neutron() ) return false;
00088
00089 if ( dic.find( std::pair < const G4Material* , const G4Element* > ( (G4Material*)NULL , element ) ) != dic.end()
00090 || dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) ) != dic.end() ) return true;
00091
00092 return false;
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102 }
00103
00104 G4double G4NeutronHPThermalScatteringData::GetIsoCrossSection( const G4DynamicParticle* dp ,
00105 G4int , G4int ,
00106 const G4Isotope* ,
00107 const G4Element* element ,
00108 const G4Material* material )
00109 {
00110 if ( dp->GetKineticEnergy() == ke_cache && element == element_cache && material == material_cache ) return xs_cache;
00111
00112 ke_cache = dp->GetKineticEnergy();
00113 element_cache = element;
00114 material_cache = material;
00115
00116 G4double xs = GetCrossSection( dp , element , material );
00117 xs_cache = xs;
00118 return xs;
00119
00120 }
00121
00122 void G4NeutronHPThermalScatteringData::clearCurrentXSData()
00123 {
00124 std::map< G4int , std::map< G4double , G4NeutronHPVector* >* >::iterator it;
00125 std::map< G4double , G4NeutronHPVector* >::iterator itt;
00126
00127 for ( it = coherent.begin() ; it != coherent.end() ; it++ )
00128 {
00129 if ( it->second != NULL )
00130 {
00131 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
00132 {
00133 delete itt->second;
00134 }
00135 }
00136 delete it->second;
00137 }
00138
00139 for ( it = incoherent.begin() ; it != incoherent.end() ; it++ )
00140 {
00141 if ( it->second != NULL )
00142 {
00143 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
00144 {
00145 delete itt->second;
00146 }
00147 }
00148 delete it->second;
00149 }
00150
00151 for ( it = inelastic.begin() ; it != inelastic.end() ; it++ )
00152 {
00153 if ( it->second != NULL )
00154 {
00155 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
00156 {
00157 delete itt->second;
00158 }
00159 }
00160 delete it->second;
00161 }
00162
00163 coherent.clear();
00164 incoherent.clear();
00165 inelastic.clear();
00166
00167 }
00168
00169
00170
00171 G4bool G4NeutronHPThermalScatteringData::IsApplicable( const G4DynamicParticle* aP , const G4Element* anEle )
00172 {
00173 G4bool result = false;
00174
00175 G4double eKin = aP->GetKineticEnergy();
00176
00177 if ( eKin < emax )
00178 {
00179
00180 if ( aP->GetDefinition() == G4Neutron::Neutron() )
00181 {
00182
00183 G4int ie = (G4int) anEle->GetIndex();
00184 std::vector < G4int >::iterator it;
00185 for ( it = indexOfThermalElement.begin() ; it != indexOfThermalElement.end() ; it++ )
00186 {
00187 if ( ie == *it ) return true;
00188 }
00189 }
00190 }
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200 return result;
00201 }
00202
00203
00204 void G4NeutronHPThermalScatteringData::BuildPhysicsTable(const G4ParticleDefinition& aP)
00205 {
00206
00207 if ( &aP != G4Neutron::Neutron() )
00208 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
00209
00210
00211 dic.clear();
00212 clearCurrentXSData();
00213 std::map < G4String , G4int > co_dic;
00214
00215
00216 static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
00217 size_t numberOfMaterials = G4Material::GetNumberOfMaterials();
00218 for ( size_t i = 0 ; i < numberOfMaterials ; i++ )
00219 {
00220 G4Material* material = (*theMaterialTable)[i];
00221 size_t numberOfElements = material->GetNumberOfElements();
00222 for ( size_t j = 0 ; j < numberOfElements ; j++ )
00223 {
00224 const G4Element* element = material->GetElement(j);
00225 if ( names->IsThisThermalElement ( material->GetName() , element->GetName() ) )
00226 {
00227 G4int ts_ID_of_this_geometry;
00228 G4String ts_ndl_name = names->GetTS_NDL_Name( material->GetName() , element->GetName() );
00229 if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
00230 {
00231 ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) -> second;
00232 }
00233 else
00234 {
00235 ts_ID_of_this_geometry = co_dic.size();
00236 co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
00237 }
00238
00239
00240
00241
00242
00243 dic.insert( std::pair < std::pair < G4Material* , const G4Element* > , G4int > ( std::pair < G4Material* , const G4Element* > ( material , element ) , ts_ID_of_this_geometry ) );
00244 }
00245 }
00246 }
00247
00248
00249 static const G4ElementTable* theElementTable = G4Element::GetElementTable();
00250 size_t numberOfElements = G4Element::GetNumberOfElements();
00251
00252 for ( size_t i = 0 ; i < numberOfElements ; i++ )
00253 {
00254 const G4Element* element = (*theElementTable)[i];
00255 if ( names->IsThisThermalElement ( element->GetName() ) )
00256 {
00257 if ( names->IsThisThermalElement ( element->GetName() ) )
00258 {
00259 G4int ts_ID_of_this_geometry;
00260 G4String ts_ndl_name = names->GetTS_NDL_Name( element->GetName() );
00261 if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
00262 {
00263 ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) -> second;
00264 }
00265 else
00266 {
00267 ts_ID_of_this_geometry = co_dic.size();
00268 co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
00269 }
00270
00271
00272
00273
00274
00275 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 ) );
00276 }
00277 }
00278 }
00279
00280 G4cout << G4endl;
00281 G4cout << "Neutron HP Thermal Scattering Data: Following material-element pairs and/or elements are registered." << G4endl;
00282 for ( std::map < std::pair < const G4Material* , const G4Element* > , G4int >::iterator it = dic.begin() ; it != dic.end() ; it++ )
00283 {
00284 if ( it->first.first != NULL )
00285 {
00286 G4cout << "Material " << it->first.first->GetName() << " - Element " << it->first.second->GetName() << ", internal thermal scattering id " << it->second << G4endl;
00287 }
00288 else
00289 {
00290 G4cout << "Element " << it->first.second->GetName() << ", internal thermal scattering id " << it->second << G4endl;
00291 }
00292 }
00293 G4cout << G4endl;
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305 G4String dirName;
00306 if ( !getenv( "G4NEUTRONHPDATA" ) )
00307 throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
00308 G4String baseName = getenv( "G4NEUTRONHPDATA" );
00309
00310 dirName = baseName + "/ThermalScattering";
00311
00312 G4String ndl_filename;
00313 G4String full_name;
00314
00315 for ( std::map < G4String , G4int >::iterator it = co_dic.begin() ; it != co_dic.end() ; it++ )
00316 {
00317 ndl_filename = it->first;
00318 G4int ts_ID = it->second;
00319
00320
00321 full_name = dirName + "/Coherent/CrossSection/" + ndl_filename;
00322 std::map< G4double , G4NeutronHPVector* >* coh_amapTemp_EnergyCross = readData( full_name );
00323 coherent.insert ( std::pair < G4int , std::map< G4double , G4NeutronHPVector* >* > ( ts_ID , coh_amapTemp_EnergyCross ) );
00324
00325
00326 full_name = dirName + "/Incoherent/CrossSection/" + ndl_filename;
00327 std::map< G4double , G4NeutronHPVector* >* incoh_amapTemp_EnergyCross = readData( full_name );
00328 incoherent.insert ( std::pair < G4int , std::map< G4double , G4NeutronHPVector* >* > ( ts_ID , incoh_amapTemp_EnergyCross ) );
00329
00330
00331 full_name = dirName + "/Inelastic/CrossSection/" + ndl_filename;
00332 std::map< G4double , G4NeutronHPVector* >* inela_amapTemp_EnergyCross = readData( full_name );
00333 inelastic.insert ( std::pair < G4int , std::map< G4double , G4NeutronHPVector* >* > ( ts_ID , inela_amapTemp_EnergyCross ) );
00334
00335 }
00336
00337 }
00338
00339
00340
00341 std::map< G4double , G4NeutronHPVector* >* G4NeutronHPThermalScatteringData::readData ( G4String full_name )
00342 {
00343
00344 std::map< G4double , G4NeutronHPVector* >* aData = new std::map< G4double , G4NeutronHPVector* >;
00345
00346 std::ifstream theChannel( full_name.c_str() );
00347
00348
00349
00350 G4int dummy;
00351 while ( theChannel >> dummy )
00352 {
00353 theChannel >> dummy;
00354 G4double temp;
00355 theChannel >> temp;
00356 G4NeutronHPVector* anEnergyCross = new G4NeutronHPVector;
00357 G4int nData;
00358 theChannel >> nData;
00359 anEnergyCross->Init ( theChannel , nData , eV , barn );
00360 aData->insert ( std::pair < G4double , G4NeutronHPVector* > ( temp , anEnergyCross ) );
00361 }
00362 theChannel.close();
00363
00364 return aData;
00365
00366 }
00367
00368
00369
00370 void G4NeutronHPThermalScatteringData::DumpPhysicsTable( const G4ParticleDefinition& aP )
00371 {
00372 if( &aP != G4Neutron::Neutron() )
00373 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
00374
00375 }
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419 G4double G4NeutronHPThermalScatteringData::GetCrossSection( const G4DynamicParticle* aP , const G4Element*anE , const G4Material* aM )
00420 {
00421 G4double result = 0;
00422
00423 G4int ts_id =getTS_ID( aM , anE );
00424
00425 if ( ts_id == -1 ) return result;
00426
00427 G4double aT = aM->GetTemperature();
00428
00429 G4double Xcoh = GetX ( aP , aT , coherent.find(ts_id)->second );
00430 G4double Xincoh = GetX ( aP , aT , incoherent.find(ts_id)->second );
00431 G4double Xinela = GetX ( aP , aT , inelastic.find(ts_id)->second );
00432
00433 result = Xcoh + Xincoh + Xinela;
00434
00435
00436
00437 return result;
00438 }
00439
00440
00441 G4double G4NeutronHPThermalScatteringData::GetInelasticCrossSection( const G4DynamicParticle* aP , const G4Element*anE , const G4Material* aM )
00442 {
00443 G4double result = 0;
00444 G4int ts_id = getTS_ID( aM , anE );
00445 G4double aT = aM->GetTemperature();
00446 result = GetX ( aP , aT , inelastic.find( ts_id )->second );
00447 return result;
00448 }
00449
00450 G4double G4NeutronHPThermalScatteringData::GetCoherentCrossSection( const G4DynamicParticle* aP , const G4Element*anE , const G4Material* aM )
00451 {
00452 G4double result = 0;
00453 G4int ts_id = getTS_ID( aM , anE );
00454 G4double aT = aM->GetTemperature();
00455 result = GetX ( aP , aT , coherent.find( ts_id )->second );
00456 return result;
00457 }
00458
00459 G4double G4NeutronHPThermalScatteringData::GetIncoherentCrossSection( const G4DynamicParticle* aP , const G4Element*anE , const G4Material* aM )
00460 {
00461 G4double result = 0;
00462 G4int ts_id = getTS_ID( aM , anE );
00463 G4double aT = aM->GetTemperature();
00464 result = GetX ( aP , aT , incoherent.find( ts_id )->second );
00465 return result;
00466 }
00467
00468
00469
00470 G4int G4NeutronHPThermalScatteringData::getTS_ID ( const G4Material* material , const G4Element* element )
00471 {
00472 G4int result = -1;
00473 if ( dic.find( std::pair < const G4Material* , const G4Element* > ( (G4Material*)NULL , element ) ) != dic.end() )
00474 return dic.find( std::pair < const G4Material* , const G4Element* > ( (G4Material*)NULL , element ) )->second;
00475 if ( dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) ) != dic.end() )
00476 return dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) )->second;
00477 return result;
00478 }
00479
00480
00481
00482
00483 G4double G4NeutronHPThermalScatteringData::GetX ( const G4DynamicParticle* aP, G4double aT , std::map < G4double , G4NeutronHPVector* >* amapTemp_EnergyCross )
00484 {
00485 G4double result = 0;
00486 if ( amapTemp_EnergyCross->size() == 0 ) return result;
00487
00488 std::map< G4double , G4NeutronHPVector* >::iterator it;
00489 for ( it = amapTemp_EnergyCross->begin() ; it != amapTemp_EnergyCross->end() ; it++ )
00490 {
00491 if ( aT < it->first ) break;
00492 }
00493 if ( it == amapTemp_EnergyCross->begin() ) it++;
00494 else if ( it == amapTemp_EnergyCross->end() ) it--;
00495
00496 G4double eKinetic = aP->GetKineticEnergy();
00497
00498 G4double TH = it->first;
00499 G4double XH = it->second->GetXsec ( eKinetic );
00500
00501
00502
00503 it--;
00504 G4double TL = it->first;
00505 G4double XL = it->second->GetXsec ( eKinetic );
00506
00507
00508
00509 if ( TH == TL )
00510 throw G4HadronicException(__FILE__, __LINE__, "Thermal Scattering Data Error!");
00511
00512 G4double T = aT;
00513 G4double X = ( XH - XL ) / ( TH - TL ) * ( T - TL ) + XL;
00514 result = X;
00515
00516 return result;
00517 }
00518
00519