G4LENDCrossSection.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 // Class Description
00027 // Cross Section for LEND (Low Energy Nuclear Data)
00028 // LEND is Geant4 interface for GIDI (General Interaction Data Interface) 
00029 // which gives a discription of nuclear and atomic reactions, such as
00030 //    Binary collision cross sections
00031 //    Particle number multiplicity distributions of reaction products
00032 //    Energy and angular distributions of reaction products
00033 //    Derived calculational constants
00034 // GIDI is developped at Lawrence Livermore National Laboratory
00035 // Class Description - End
00036 
00037 // 071025 First implementation done by T. Koi (SLAC/SCCS)
00038 // 101118 Name modifications for release T. Koi (SLAC/PPA)
00039 
00040 #include "G4LENDCrossSection.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "G4ElementTable.hh"
00043 #include "G4HadronicException.hh"
00044 
00045 //TK110811
00046 G4bool G4LENDCrossSection::IsIsoApplicable( const G4DynamicParticle* dp, G4int /*iZ*/ , G4int /*aA*/ , 
00047                                             const G4Element* /*element*/ , const G4Material* /*material*/ )
00048 {
00049    G4double eKin = dp->GetKineticEnergy();
00050    if ( dp->GetDefinition() != proj ) return false;
00051    if ( eKin > GetMaxKinEnergy() || eKin < GetMinKinEnergy() ) return false;
00052 
00053    return true;
00054 }
00055 
00056 G4double G4LENDCrossSection::GetIsoCrossSection( const G4DynamicParticle* dp , G4int iZ , G4int iA ,
00057                                                  const G4Isotope* /*isotope*/ , const G4Element* /*elment*/ , const G4Material* material )
00058 {
00059 
00060    G4double xs = 0.0;
00061    G4double ke = dp->GetKineticEnergy();
00062    G4double temp = material->GetTemperature();
00063 
00064    G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA ) )->second->GetTarget();
00065 
00066    xs = getLENDCrossSection ( aTarget , ke , temp );
00067 
00068    return xs;
00069 }
00070 
00071 
00072 /*
00073 G4bool G4LENDCrossSection::IsApplicable(const G4DynamicParticle*aP, const G4Element*)
00074 {
00075    G4bool result = true;
00076    G4double eKin = aP->GetKineticEnergy();
00077    if( eKin > GetMaxKinEnergy() || aP->GetDefinition() != proj ) result = false;
00078    return result;
00079 }
00080 */
00081 
00082 G4LENDCrossSection::G4LENDCrossSection( const G4String nam )
00083 :G4VCrossSectionDataSet( nam )
00084 {
00085 
00086    default_evaluation = "endl99";
00087    allow_nat = false;
00088    allow_any = false;
00089 
00090    SetMinKinEnergy(  0*MeV );
00091    SetMaxKinEnergy( 20*MeV );
00092 
00093    lend_manager = G4LENDManager::GetInstance(); 
00094 
00095 }
00096    
00097 G4LENDCrossSection::~G4LENDCrossSection()
00098 {
00099 
00100    for ( std::map< G4int , G4LENDUsedTarget* >::iterator 
00101          it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
00102    { 
00103       delete it->second;  
00104    }
00105 
00106 }
00107    
00108 void G4LENDCrossSection::BuildPhysicsTable( const G4ParticleDefinition&  )
00109 {
00110    create_used_target_map();
00111 }
00112 
00113 void G4LENDCrossSection::DumpPhysicsTable(const G4ParticleDefinition& aP)
00114 {
00115 
00116   if ( &aP != proj ) 
00117      throw G4HadronicException(__FILE__, __LINE__, "Attempt to use LEND data for particles other than neutrons!!!");  
00118 
00119    G4cout << G4endl;
00120    G4cout << "Dump Cross Sections of " << GetName() << G4endl;
00121    G4cout << "(Pointwise cross-section at 300 Kelvin.)" << G4endl;
00122    G4cout << G4endl;
00123 
00124    G4cout << "Target informaiton " << G4endl;
00125 
00126    for ( std::map< G4int , G4LENDUsedTarget* >::iterator 
00127          it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
00128    {
00129       G4cout 
00130          << "Wanted " << it->second->GetWantedEvaluation() 
00131          << ", Z= " << it->second->GetWantedZ() 
00132          << ", A= " << it->second->GetWantedA() 
00133          << "; Actual " << it->second->GetActualEvaluation() 
00134          << ", Z= " << it->second->GetActualZ() 
00135          << ", A= " << it->second->GetActualA() 
00136          << ", " << it->second->GetTarget() 
00137          << G4endl; 
00138 
00139       G4int ie = 0;
00140 
00141       G4GIDI_target* aTarget = it->second->GetTarget();
00142       G4double aT = 300;
00143       for ( ie = 0 ; ie < 130 ; ie++ )
00144       {
00145          G4double ke = 1.0e-5 * std::pow ( 10.0 , ie/10.0 ) *eV;
00146 
00147          if ( ke < 20*MeV )
00148          {
00149             G4cout << "  " << GetName() << ", cross section at " << ke/eV << " [eV] = " << getLENDCrossSection ( aTarget , ke , aT )/barn << " [barn] " << G4endl;
00150          }
00151       }
00152       G4cout << G4endl;
00153 
00154    }
00155 
00156 }
00157 
00158 
00159 /*
00160 //110810
00161 //G4double G4LENDCrossSection::GetCrossSection(const G4DynamicParticle* aP , const G4Element* anElement , G4double aT)
00162 G4double G4LENDCrossSection::GetCrossSection(const G4DynamicParticle* aP , int iZ , const G4Material* aMat)
00163 {
00164 
00165 //110810
00166    G4double aT = aMat->GetTemperature();
00167    G4Element* anElement = lend_manager->GetNistElementBuilder()->FindOrBuildElement( iZ );
00168 
00169    G4double ke = aP->GetKineticEnergy();
00170    G4double XS = 0.0;
00171 
00172    G4int numberOfIsotope = anElement->GetNumberOfIsotopes(); 
00173 
00174    if ( numberOfIsotope > 0 )
00175    {
00176       // User Defined Abundances   
00177       for ( G4int i_iso = 0 ; i_iso < numberOfIsotope ; i_iso++ )
00178       {
00179 
00180          G4int iZ = anElement->GetIsotope( i_iso )->GetZ();
00181          G4int iA = anElement->GetIsotope( i_iso )->GetN();
00182          G4double ratio = anElement->GetRelativeAbundanceVector()[i_iso];
00183 
00184          G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA ) )->second->GetTarget();
00185          XS += ratio*getLENDCrossSection ( aTarget , ke , aT );
00186 
00187       }
00188    }
00189    else
00190    {
00191       // Natural Abundances   
00192       G4NistElementBuilder* nistElementBuild = lend_manager->GetNistElementBuilder();
00193       G4int iZ = int ( anElement->GetZ() );
00194       G4int numberOfNistIso = nistElementBuild->GetNumberOfNistIsotopes( int ( anElement->GetZ() ) ); 
00195 
00196        G4int Nfirst = nistElementBuild->GetNistFirstIsotopeN( iZ );
00197       for ( G4int i = 0 ; i < numberOfNistIso ; i++ )
00198       {
00199          G4int iA = Nfirst + i;  
00200          G4double ratio = nistElementBuild->GetIsotopeAbundance( iZ , iA );
00201          if ( ratio > 0.0 )
00202          {
00203             G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA ) )->second->GetTarget();
00204             XS += ratio*getLENDCrossSection ( aTarget , ke , aT );
00205             //G4cout << ke/eV << " "  << iZ << " " << iMass << " " << aTarget << " " << getLENDCrossSection ( aTarget , ke , aT ) << G4endl;
00206          }
00207       }
00208    }
00209  
00210    //G4cout << "XS= " << XS << G4endl;
00211    return XS;
00212 }
00213 
00214 
00215 
00216 //110810
00217 //G4double G4LENDCrossSection::GetIsoCrossSection(const G4DynamicParticle* dp, const G4Isotope* isotope, G4double aT )
00218 G4double G4LENDCrossSection::GetIsoCrossSection(const G4DynamicParticle* dp, const G4Isotope* isotope, const G4Material* aMat)
00219 {
00220 
00221 //110810
00222    G4double aT = aMat->GetTemperature();
00223 
00224    G4double ke = dp->GetKineticEnergy();
00225 
00226    G4int iZ = isotope->GetZ();
00227    G4int iA = isotope->GetN();
00228 
00229    G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA ) )->second->GetTarget();
00230 
00231    return getLENDCrossSection ( aTarget , ke , aT );
00232 
00233 }
00234 
00235 
00236 
00237 //110810
00238 //G4double G4LENDCrossSection::GetZandACrossSection(const G4DynamicParticle* dp, G4int iZ, G4int iA, G4double aT)
00239 G4double G4LENDCrossSection::GetZandACrossSection(const G4DynamicParticle* dp, G4int iZ, G4int iA, const G4Material* aMat)
00240 {
00241 
00242 //110810
00243    G4double aT = aMat->GetTemperature();
00244 
00245    G4double ke = dp->GetKineticEnergy();
00246 
00247    G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA ) )->second->GetTarget();
00248 
00249    return getLENDCrossSection ( aTarget , ke , aT );
00250 
00251 }
00252 */
00253 
00254 
00255 
00256 void G4LENDCrossSection::recreate_used_target_map()
00257 {
00258    for ( std::map< G4int , G4LENDUsedTarget* >::iterator 
00259          it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
00260    { 
00261       delete it->second;  
00262    }
00263    usedTarget_map.clear();
00264 
00265    create_used_target_map();
00266 }
00267 
00268 
00269 
00270 void G4LENDCrossSection::create_used_target_map()
00271 {
00272 
00273    lend_manager->RequestChangeOfVerboseLevel( verboseLevel );
00274 
00275    size_t numberOfElements = G4Element::GetNumberOfElements();
00276    static const G4ElementTable* theElementTable = G4Element::GetElementTable();
00277 
00278    for ( size_t i = 0 ; i < numberOfElements ; ++i )
00279    {
00280 
00281       const G4Element* anElement = (*theElementTable)[i];
00282       G4int numberOfIsotope = anElement->GetNumberOfIsotopes(); 
00283 
00284       if ( numberOfIsotope > 0 )
00285       {
00286       // User Defined Abundances   
00287          for ( G4int i_iso = 0 ; i_iso < numberOfIsotope ; i_iso++ )
00288          {
00289             G4int iZ = anElement->GetIsotope( i_iso )->GetZ();
00290             G4int iA = anElement->GetIsotope( i_iso )->GetN();
00291 
00292             //G4LENDUsedTarget* aTarget = new G4LENDUsedTarget ( G4Neutron::Neutron() , default_evaluation , iZ , iA );  
00293             G4LENDUsedTarget* aTarget = new G4LENDUsedTarget ( proj , default_evaluation , iZ , iA );  
00294             if ( allow_nat == true ) aTarget->AllowNat();
00295             if ( allow_any == true ) aTarget->AllowAny();
00296             usedTarget_map.insert( std::pair< G4int , G4LENDUsedTarget* > ( lend_manager->GetNucleusEncoding( iZ , iA ) , aTarget ) );
00297          }
00298       }
00299       else
00300       {
00301       // Natural Abundances   
00302          G4NistElementBuilder* nistElementBuild = lend_manager->GetNistElementBuilder();
00303          G4int iZ = int ( anElement->GetZ() );
00304          //G4cout << nistElementBuild->GetNumberOfNistIsotopes( int ( anElement->GetZ() ) ) << G4endl;
00305          G4int numberOfNistIso = nistElementBuild->GetNumberOfNistIsotopes( int ( anElement->GetZ() ) ); 
00306 
00307          for ( G4int ii = 0 ; ii < numberOfNistIso ; ii++ )
00308          {
00309             //G4cout << nistElementBuild->GetIsotopeAbundance( iZ , nistElementBuild->GetNistFirstIsotopeN( iZ ) + i ) << G4endl;
00310             if ( nistElementBuild->GetIsotopeAbundance( iZ , nistElementBuild->GetNistFirstIsotopeN( iZ ) + ii ) > 0 )
00311             {
00312                G4int iMass = nistElementBuild->GetNistFirstIsotopeN( iZ ) + ii;  
00313                //G4cout << iZ << " " << nistElementBuild->GetNistFirstIsotopeN( iZ ) + i << " " << nistElementBuild->GetIsotopeAbundance ( iZ , iMass ) << G4endl;  
00314 
00315                G4LENDUsedTarget* aTarget = new G4LENDUsedTarget ( proj , default_evaluation , iZ , iMass );  
00316                if ( allow_nat == true ) aTarget->AllowNat();
00317                if ( allow_any == true ) aTarget->AllowAny();
00318                usedTarget_map.insert( std::pair< G4int , G4LENDUsedTarget* > ( lend_manager->GetNucleusEncoding( iZ , iMass ) , aTarget ) );
00319 
00320             }
00321 
00322          }
00323       }
00324    }
00325 
00326    G4cout << "Dump UsedTarget for " << GetName() << G4endl;
00327    G4cout << "Requested Evaluation, Z , A -> Actual Evaluation, Z , A(0=Nat) , Pointer of Target" << G4endl;
00328    for ( std::map< G4int , G4LENDUsedTarget* >::iterator 
00329          it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
00330    {
00331       G4cout 
00332          << " " << it->second->GetWantedEvaluation() 
00333          << ", " << it->second->GetWantedZ() 
00334          << ", " << it->second->GetWantedA() 
00335          << " -> " << it->second->GetActualEvaluation() 
00336          << ", " << it->second->GetActualZ() 
00337          << ", " << it->second->GetActualA() 
00338          << ", " << it->second->GetTarget() 
00339          << G4endl; 
00340    } 
00341 
00342 }
00343 
00344                                                            // elow          ehigh       xs_elow      xs_ehigh      ke (ke < elow)
00345 G4double G4LENDCrossSection::GetUltraLowEnergyExtrapolatedXS( G4double x1, G4double x2, G4double y1, G4double y2 , G4double ke )
00346 {
00347    //XS propotinal to 1/v at low energy -> 1/root(E) 
00348    //XS = a * 1/root(E) + b  
00349    G4double a = ( y2 - y1 ) / ( 1/std::sqrt(x2) - 1/std::sqrt(x1) );
00350    G4double b = y1 - a * 1/std::sqrt(x1);
00351    G4double result = a * 1/std::sqrt(ke) + b;
00352    return result;
00353 }

Generated on Mon May 27 17:48:45 2013 for Geant4 by  doxygen 1.4.7