G4NeutronHPThermalScatteringData.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // Thermal Neutron Scattering
00027 // Koi, Tatsumi (SCCS/SLAC)
00028 //
00029 // Class Description
00030 // Cross Sections for a high precision (based on evaluated data
00031 // libraries) description of themal neutron scattering below 4 eV;
00032 // Based on Thermal neutron scattering files
00033 // from the evaluated nuclear data files ENDF/B-VI, Release2
00034 // To be used in your physics list in case you need this physics.
00035 // In this case you want to register an object of this class with
00036 // the corresponding process.
00037 // Class Description - End
00038 
00039 // 15-Nov-06 First implementation is done by T. Koi (SLAC/SCCS)
00040 // 070625 implement clearCurrentXSData to fix memory leaking by T. Koi
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 //#include "G4NeutronHPData.hh"
00050 
00051 G4NeutronHPThermalScatteringData::G4NeutronHPThermalScatteringData()
00052 :G4VCrossSectionDataSet("NeutronHPThermalScatteringData")
00053 {
00054 // Upper limit of neutron energy 
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    //BuildPhysicsTable( *G4Neutron::Neutron() );
00069 }
00070 
00071 G4NeutronHPThermalScatteringData::~G4NeutronHPThermalScatteringData()
00072 {
00073 
00074    clearCurrentXSData();
00075 
00076    delete names;
00077 }
00078 
00079 G4bool G4NeutronHPThermalScatteringData::IsIsoApplicable( const G4DynamicParticle* dp , 
00080                                                 G4int /*Z*/ , G4int /*A*/ ,
00081                                                 const G4Element* element ,
00082                                                 const G4Material* material )
00083 {
00084    G4double eKin = dp->GetKineticEnergy();
00085    if ( eKin > 4.0*eV //GetMaxKinEnergy() 
00086      || eKin < 0 //GetMinKinEnergy() 
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 //   return IsApplicable( dp , element );
00095 /*
00096    G4double eKin = dp->GetKineticEnergy();
00097    if ( eKin > 4.0*eV //GetMaxKinEnergy() 
00098      || eKin < 0 //GetMinKinEnergy() 
00099      || dp->GetDefinition() != G4Neutron::Neutron() ) return false;                                   
00100    return true;
00101 */
00102 }
00103 
00104 G4double G4NeutronHPThermalScatteringData::GetIsoCrossSection( const G4DynamicParticle* dp ,
00105                                    G4int /*Z*/ , G4int /*A*/ ,
00106                                    const G4Isotope* /*iso*/  ,
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    //G4double xs = GetCrossSection( dp , element , material->GetTemperature() );
00116    G4double xs = GetCrossSection( dp , element , material );
00117    xs_cache = xs;
00118    return xs;
00119    //return GetCrossSection( dp , element , material->GetTemperature() );
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    // Check energy 
00177    if ( eKin < emax )
00178    {
00179       // Check Particle Species
00180       if ( aP->GetDefinition() == G4Neutron::Neutron() ) 
00181       {
00182         // anEle is one of Thermal elements 
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    if ( names->IsThisThermalElement ( anEle->GetName() ) )
00194    {
00195       // Check energy and projectile species 
00196       G4double eKin = aP->GetKineticEnergy();
00197       if ( eKin < emax && aP->GetDefinition() == G4Neutron::Neutron() ) result = true; 
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    //std::map < std::pair < G4Material* , const G4Element* > , G4int > dic;   
00211    dic.clear();   
00212    clearCurrentXSData();
00213    std::map < G4String , G4int > co_dic;   
00214 
00215    //Searching Nist Materials
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             //G4cout << "Neutron HP Thermal Scattering Data : Registering a material-element pair of " 
00240             //       << material->GetName() << " " << element->GetName() 
00241             //       << " as internal thermal scattering id of  " <<  ts_ID_of_this_geometry << "." << G4endl;
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    //Searching TS Elements 
00249    static const G4ElementTable* theElementTable = G4Element::GetElementTable();
00250    size_t numberOfElements = G4Element::GetNumberOfElements();
00251    //size_t numberOfThermalElements = 0; 
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             //G4cout << "Neutron HP Thermal Scattering: Registering an element of " 
00272             //       << material->GetName() << " " << element->GetName() 
00273             //       << " as internal thermal scattering id of  " <<  ts_ID_of_this_geometry << "." << G4endl;
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    //G4cout << "Neutron HP Thermal Scattering Data: Following NDL thermal scattering files are assigned to the internal thermal scattering id." << G4endl;
00297    //for ( std::map < G4String , G4int >::iterator it = co_dic.begin() ; it != co_dic.end() ; it++ )  
00298    //{
00299    //   G4cout << "NDL file name " << it->first << ", internal thermal scattering id " << it->second << G4endl;
00300    //}
00301 
00302 
00303    // Read Cross Section Data files
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       // Coherent
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       // Incoherent
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       // Inelastic
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    //G4cout << "G4NeutronHPThermalScatteringData " << name << G4endl;
00349 
00350    G4int dummy; 
00351    while ( theChannel >> dummy )   // MF
00352    {
00353       theChannel >> dummy;   // MT
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 //  G4cout << "G4NeutronHPThermalScatteringData::DumpPhysicsTable still to be implemented"<<G4endl;
00375 }
00376 
00377 //#include "G4Nucleus.hh"
00378 //#include "G4NucleiPropertiesTable.hh"
00379 //#include "G4Neutron.hh"
00380 //#include "G4Electron.hh"
00381 
00382 
00383 
00384 /*
00385 G4double G4NeutronHPThermalScatteringData::GetCrossSection( const G4DynamicParticle* aP , const G4Element*anE , G4double aT )
00386 {
00387 
00388    G4double result = 0;
00389    const G4Material* aM = NULL;
00390 
00391    G4int iele = anE->GetIndex();
00392 
00393    if ( dic.find( std::pair < const G4Material* , const G4Element* > ( (G4Material*)NULL , anE ) ) != dic.end() )
00394    {
00395       iele = dic.find( std::pair < const G4Material* , const G4Element* > ( (G4Material*)NULL , anE ) )->second;
00396    }
00397    else if ( dic.find( std::pair < const G4Material* , const G4Element* > ( aM , anE ) ) != dic.end() )
00398    {
00399       iele = dic.find( std::pair < const G4Material* , const G4Element* > ( aM , anE ) )->second;
00400    }
00401    else
00402    {
00403       return result;
00404    }
00405 
00406    G4double Xcoh = GetX ( aP , aT , coherent.find(iele)->second );
00407    G4double Xincoh = GetX ( aP , aT , incoherent.find(iele)->second );
00408    G4double Xinela = GetX ( aP , aT , inelastic.find(iele)->second );
00409 
00410    result = Xcoh + Xincoh + Xinela;
00411 
00412    //G4cout << "G4NeutronHPThermalScatteringData::GetCrossSection  Tot= " << result/barn << " Coherent= " << Xcoh/barn << " Incoherent= " << Xincoh/barn << " Inelastic= " << Xinela/barn << G4endl;
00413 
00414    return result;
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    //G4cout << "G4NeutronHPThermalScatteringData::GetCrossSection  Tot= " << result/barn << " Coherent= " << Xcoh/barn << " Incoherent= " << Xincoh/barn << " Inelastic= " << Xinela/barn << G4endl;
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++;  // lower than first
00494       else if ( it == amapTemp_EnergyCross->end() ) it--;  // upper than last
00495 
00496    G4double eKinetic = aP->GetKineticEnergy();
00497 
00498    G4double TH = it->first;
00499    G4double XH = it->second->GetXsec ( eKinetic ); 
00500 
00501    //G4cout << "G4NeutronHPThermalScatteringData::GetX TH " << TH << " E " << eKinetic <<  " XH " << XH << G4endl;
00502 
00503    it--;
00504    G4double TL = it->first;
00505    G4double XL = it->second->GetXsec ( eKinetic ); 
00506 
00507    //G4cout << "G4NeutronHPThermalScatteringData::GetX TL " << TL << " E " << eKinetic <<  " XL " << XL << G4endl;
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 

Generated on Mon May 27 17:49:04 2013 for Geant4 by  doxygen 1.4.7