G4NeutronHPCaptureData.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 // neutron_hp -- source file
00027 // J.P. Wellisch, Nov-1996
00028 // A prototype of the low energy neutron transport model.
00029 //
00030 // 070523 add neglecting doppler broadening on the fly. T. Koi
00031 // 070613 fix memory leaking by T. Koi
00032 // 071002 enable cross section dump by T. Koi
00033 // 080428 change checking point of "neglecting doppler broadening" flag 
00034 //        from GetCrossSection to BuildPhysicsTable by T. Koi
00035 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
00036 //
00037 #include "G4NeutronHPCaptureData.hh"
00038 #include "G4PhysicalConstants.hh"
00039 #include "G4SystemOfUnits.hh"
00040 #include "G4Neutron.hh"
00041 #include "G4ElementTable.hh"
00042 #include "G4NeutronHPData.hh"
00043 
00044 G4NeutronHPCaptureData::G4NeutronHPCaptureData()
00045 :G4VCrossSectionDataSet("NeutronHPCaptureXS")
00046 {
00047    SetMinKinEnergy( 0*MeV );                                   
00048    SetMaxKinEnergy( 20*MeV );                                   
00049 
00050    ke_cache = 0.0;
00051    xs_cache = 0.0;
00052    element_cache = NULL;
00053    material_cache = NULL;
00054 
00055    theCrossSections = 0;
00056    onFlightDB = true;
00057 
00058    BuildPhysicsTable(*G4Neutron::Neutron());
00059 }
00060    
00061 G4NeutronHPCaptureData::~G4NeutronHPCaptureData()
00062 {
00063    if ( theCrossSections != 0 ) theCrossSections->clearAndDestroy();
00064 
00065    delete theCrossSections;
00066 }
00067    
00068 G4bool G4NeutronHPCaptureData::IsIsoApplicable( const G4DynamicParticle* dp , 
00069                                                 G4int /*Z*/ , G4int /*A*/ ,
00070                                                 const G4Element* /*elm*/ ,
00071                                                 const G4Material* /*mat*/ )
00072 {
00073    G4double eKin = dp->GetKineticEnergy();
00074    if ( eKin > GetMaxKinEnergy() 
00075      || eKin < GetMinKinEnergy() 
00076      || dp->GetDefinition() != G4Neutron::Neutron() ) return false;                                   
00077 
00078    return true;
00079 }
00080 
00081 G4double G4NeutronHPCaptureData::GetIsoCrossSection( const G4DynamicParticle* dp ,
00082                                    G4int /*Z*/ , G4int /*A*/ ,
00083                                    const G4Isotope* /*iso*/  ,
00084                                    const G4Element* element ,
00085                                    const G4Material* material )
00086 {
00087    if ( dp->GetKineticEnergy() == ke_cache && element == element_cache &&  material == material_cache ) return xs_cache;
00088 
00089    ke_cache = dp->GetKineticEnergy();
00090    element_cache = element;
00091    material_cache = material;
00092    G4double xs = GetCrossSection( dp , element , material->GetTemperature() );
00093    xs_cache = xs;
00094    return xs;
00095 }
00096 
00097 /*
00098 G4bool G4NeutronHPCaptureData::IsApplicable(const G4DynamicParticle*aP, const G4Element*)
00099 {
00100   G4bool result = true;
00101   G4double eKin = aP->GetKineticEnergy();
00102   if(eKin>20*MeV||aP->GetDefinition()!=G4Neutron::Neutron()) result = false;
00103   return result;
00104 }
00105 */
00106 
00107 void G4NeutronHPCaptureData::BuildPhysicsTable(const G4ParticleDefinition& aP)
00108 {
00109   if(&aP!=G4Neutron::Neutron()) 
00110      throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");  
00111 
00112 //080428
00113    if ( getenv( "G4NEUTRONHP_NEGLECT_DOPPLER" ) ) 
00114    {
00115       G4cout << "Find environment variable of \"G4NEUTRONHP_NEGLECT_DOPPLER\"." << G4endl;
00116       G4cout << "On the fly Doppler broadening will be neglect in the cross section calculation of capture reaction of neutrons (<20MeV)." << G4endl;
00117       onFlightDB = false;
00118    }
00119   
00120   size_t numberOfElements = G4Element::GetNumberOfElements();
00121   // G4cout << "CALLED G4NeutronHPCaptureData::BuildPhysicsTable "<<numberOfElements<<G4endl;
00122    // TKDB
00123    //if ( theCrossSections == 0 ) theCrossSections = new G4PhysicsTable( numberOfElements );
00124    if ( theCrossSections == NULL ) 
00125       theCrossSections = new G4PhysicsTable( numberOfElements );
00126    else
00127       theCrossSections->clearAndDestroy();
00128 
00129   // make a PhysicsVector for each element
00130 
00131   static const G4ElementTable *theElementTable = G4Element::GetElementTable();
00132   for( size_t i=0; i<numberOfElements; ++i )
00133   {
00134      if(getenv("CaptureDataIndexDebug"))
00135      {
00136        G4int index_debug = ((*theElementTable)[i])->GetIndex();
00137        G4cout << "IndexDebug "<< i <<" "<<index_debug<<G4endl;
00138      }
00139      G4PhysicsVector* physVec = G4NeutronHPData::
00140       Instance()->MakePhysicsVector((*theElementTable)[i], this);
00141      theCrossSections->push_back(physVec);
00142   }
00143 }
00144 
00145 void G4NeutronHPCaptureData::DumpPhysicsTable(const G4ParticleDefinition& aP)
00146 {
00147   if(&aP!=G4Neutron::Neutron()) 
00148      throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");  
00149 
00150 //
00151 // Dump element based cross section
00152 // range 10e-5 eV to 20 MeV
00153 // 10 point per decade
00154 // in barn
00155 //
00156 
00157    G4cout << G4endl;
00158    G4cout << G4endl;
00159    G4cout << "Capture Cross Section of Neutron HP"<< G4endl;
00160    G4cout << "(Pointwise cross-section at 0 Kelvin.)" << G4endl;
00161    G4cout << G4endl;
00162    G4cout << "Name of Element" << G4endl;
00163    G4cout << "Energy[eV]  XS[barn]" << G4endl;
00164    G4cout << G4endl;
00165 
00166    size_t numberOfElements = G4Element::GetNumberOfElements();
00167    static const G4ElementTable *theElementTable = G4Element::GetElementTable();
00168 
00169    for ( size_t i = 0 ; i < numberOfElements ; ++i )
00170    {
00171 
00172       G4cout << (*theElementTable)[i]->GetName() << G4endl;
00173 
00174       G4int ie = 0;
00175 
00176       for ( ie = 0 ; ie < 130 ; ie++ )
00177       {
00178          G4double eKinetic = 1.0e-5 * std::pow ( 10.0 , ie/10.0 ) *eV;
00179          G4bool outOfRange = false;
00180 
00181          if ( eKinetic < 20*MeV )
00182          {
00183             G4cout << eKinetic/eV << " " << (*((*theCrossSections)(i))).GetValue(eKinetic, outOfRange)/barn << G4endl;
00184          }
00185 
00186       }
00187 
00188       G4cout << G4endl;
00189    }
00190 
00191 
00192 //  G4cout << "G4NeutronHPCaptureData::DumpPhysicsTable still to be implemented"<<G4endl;
00193 }
00194 
00195 #include "G4NucleiProperties.hh"
00196 
00197 G4double G4NeutronHPCaptureData::
00198 GetCrossSection(const G4DynamicParticle* aP, const G4Element*anE, G4double aT)
00199 {
00200   G4double result = 0;
00201   G4bool outOfRange;
00202   G4int index = anE->GetIndex();
00203 
00204   // prepare neutron
00205   G4double eKinetic = aP->GetKineticEnergy();
00206 
00207 //if ( getenv( "G4NEUTRONHP_NEGLECT_DOPPLER" ) )
00208 //080428
00209   if ( !onFlightDB )
00210   {
00211      G4double factor = 1.0;
00212      if ( eKinetic < aT * k_Boltzmann ) 
00213      {
00214         // below 0.1 eV neutrons 
00215         // Have to do some, but now just igonre.   
00216         // Will take care after performance check.  
00217         // factor = factor * targetV;
00218      }
00219      return ( (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) )* factor; 
00220   }
00221 
00222   G4ReactionProduct theNeutron( aP->GetDefinition() );
00223   theNeutron.SetMomentum( aP->GetMomentum() );
00224   theNeutron.SetKineticEnergy( eKinetic );
00225 
00226   // prepare thermal nucleus
00227   G4Nucleus aNuc;
00228   G4double eps = 0.0001;
00229   G4double theA = anE->GetN();
00230   G4double theZ = anE->GetZ();
00231   G4double eleMass; 
00232   eleMass = G4NucleiProperties::GetNuclearMass( static_cast<G4int>(theA+eps) , static_cast<G4int>(theZ+eps) ) / G4Neutron::Neutron()->GetPDGMass();
00233   
00234   G4ReactionProduct boosted;
00235   G4double aXsection;
00236   
00237   // MC integration loop
00238   G4int counter = 0;
00239   G4double buffer = 0;
00240   G4int size = G4int(std::max(10., aT/60*kelvin));
00241   G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutron.GetMomentum();
00242   G4double neutronVMag = neutronVelocity.mag();
00243 
00244   while(counter == 0 || std::abs(buffer-result/std::max(1,counter)) > 0.03*buffer)
00245   {
00246     if(counter) buffer = result/counter;
00247     while (counter<size)
00248     {
00249       counter ++;
00250       G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus(eleMass, aT);
00251       boosted.Lorentz(theNeutron, aThermalNuc);
00252       G4double theEkin = boosted.GetKineticEnergy();
00253       aXsection = (*((*theCrossSections)(index))).GetValue(theEkin, outOfRange);
00254       // velocity correction, or luminosity factor...
00255       G4ThreeVector targetVelocity = 1./aThermalNuc.GetMass()*aThermalNuc.GetMomentum();
00256       aXsection *= (targetVelocity-neutronVelocity).mag()/neutronVMag;
00257       result += aXsection;
00258     }
00259     size += size;
00260   }
00261   result /= counter;
00262 /*
00263   // Checking impact of  G4NEUTRONHP_NEGLECT_DOPPLER
00264   G4cout << " result " << result << " " 
00265          << (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) << " " 
00266          << (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) /result << G4endl;
00267 */
00268   return result;
00269 }

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