G4NeutronHPorLEInelasticData.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 //
00027 // 05-11-21 NeutronHP or Low Energy Parameterization Models 
00028 //          Implemented by T. Koi (SLAC/SCCS)
00029 //          If NeutronHP data do not available for an element, then Low Energy 
00030 //          Parameterization models handle the interactions of the element.
00031 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
00032 //
00033 
00034 #include "G4NeutronHPorLEInelasticData.hh"
00035 #include "G4SystemOfUnits.hh"
00036 #include "G4Neutron.hh"
00037 #include "G4ElementTable.hh"
00038 #include "G4NeutronHPData.hh"
00039 
00040 #include "G4PhysicsVector.hh"
00041 
00042 G4NeutronHPorLEInelasticData::G4NeutronHPorLEInelasticData()
00043 {
00044    SetMinKinEnergy( 0*MeV );                                   
00045    SetMaxKinEnergy( 20*MeV );                                   
00046 
00047    ke_cache = 0.0;
00048    xs_cache = 0.0;
00049    element_cache = NULL;
00050    material_cache = NULL;
00051 //   BuildPhysicsTable(*G4Neutron::Neutron());
00052 }
00053 
00054 G4NeutronHPorLEInelasticData::~G4NeutronHPorLEInelasticData()
00055 {
00056 //  delete theCrossSections;
00057 }
00058 
00059 G4bool G4NeutronHPorLEInelasticData::IsIsoApplicable( const G4DynamicParticle* dp , 
00060                                                 G4int /*Z*/ , G4int /*A*/ ,
00061                                                 const G4Element* element ,
00062                                                 const G4Material* /*mat*/ )
00063 {
00064    G4double eKin = dp->GetKineticEnergy();
00065    if ( eKin > GetMaxKinEnergy() 
00066      || eKin < GetMinKinEnergy() 
00067      || dp->GetDefinition() != G4Neutron::Neutron() ) return false;                                   
00068    if ( unavailable_elements->find( element->GetName() ) != unavailable_elements->end() ) return false;
00069 
00070    return true;
00071 }
00072 
00073 G4double G4NeutronHPorLEInelasticData::GetIsoCrossSection( const G4DynamicParticle* dp ,
00074                                    G4int /*Z*/ , G4int /*A*/ ,
00075                                    const G4Isotope* /*iso*/  ,
00076                                    const G4Element* element ,
00077                                    const G4Material* material )
00078 {
00079    if ( dp->GetKineticEnergy() == ke_cache && element == element_cache &&  material == material_cache ) return xs_cache;
00080 
00081    ke_cache = dp->GetKineticEnergy();
00082    element_cache = element;
00083    material_cache = material;
00084    G4double xs = GetCrossSection( dp , element , material->GetTemperature() );
00085    xs_cache = xs;
00086    return xs;
00087    //return GetCrossSection( dp , element , material->GetTemperature() );
00088 }
00089 
00090 G4NeutronHPorLEInelasticData::G4NeutronHPorLEInelasticData( G4NeutronHPChannelList* pChannel , std::set< G4String >* pSet )
00091 :G4VCrossSectionDataSet("NeutronHPorLEInelasticXS")
00092 {
00093    theInelasticChannel = pChannel;
00094    unavailable_elements = pSet;   
00095    //BuildPhysicsTable(*G4Neutron::Neutron());
00096 
00097    SetMinKinEnergy( 0*MeV );                                   
00098    SetMaxKinEnergy( 20*MeV );                                   
00099 
00100    ke_cache = 0.0;
00101    xs_cache = 0.0;
00102    element_cache = NULL;
00103    material_cache = NULL;
00104 }
00105  
00106 /*
00107 G4bool G4NeutronHPorLEInelasticData::IsApplicable(const G4DynamicParticle*aP, const G4Element* anElement)
00108 {
00109    G4bool result = true;
00110    G4double eKin = aP->GetKineticEnergy();
00111    if(eKin>20*MeV||aP->GetDefinition()!=G4Neutron::Neutron()) result = false;
00112    if ( unavailable_elements->find( anElement->GetName() ) != unavailable_elements->end() ) result = false;
00113    return result;
00114 }
00115 */
00116 
00117 #include "G4NeutronHPInelasticData.hh"
00118 #include "G4LPhysicsFreeVector.hh"
00119 //#include "G4NeutronHPElementData.hh"
00120    
00121 void G4NeutronHPorLEInelasticData::BuildPhysicsTable( const G4ParticleDefinition& aP )
00122 {
00123    if( &aP!=G4Neutron::Neutron() ) 
00124       throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");  
00125 
00126    size_t numberOfElements = G4Element::GetNumberOfElements();
00127    theCrossSections = new G4PhysicsTable( numberOfElements );
00128 
00129    static const G4ElementTable *theElementTable = G4Element::GetElementTable(); 
00130    for ( size_t i=0 ; i < numberOfElements; ++i )
00131    {
00132       G4PhysicsVector* thePhysVec = new G4LPhysicsFreeVector(0, 0, 0);
00133 
00134       if ( unavailable_elements->find( (*theElementTable)[i]->GetName() ) == unavailable_elements->end() ) 
00135       { 
00136 
00137          G4NeutronHPElementData* theElementData = new G4NeutronHPElementData();
00138          theElementData->Init( (*theElementTable)[i] );  
00139 
00140          G4NeutronHPVector* theHPVector = theElementData->GetData( (G4NeutronHPInelasticData*)this ); 
00141 
00142          G4int len = theHPVector->GetVectorLength();
00143 
00144          if ( len!=0 ) 
00145          {
00146             G4double emin = theHPVector->GetX(0);
00147             G4double emax = theHPVector->GetX(len-1);
00148 
00149             G4LPhysicsFreeVector* aPhysVector= new G4LPhysicsFreeVector ( len , emin , emax );
00150             for ( G4int ii=0; ii<len; ii++ )
00151             {
00152                aPhysVector->PutValues( ii , theHPVector->GetX(ii) , theHPVector->GetY(ii) );
00153             }
00154             delete thePhysVec;
00155             thePhysVec = aPhysVector;
00156          }
00157 
00158          //G4PhysicsVector* physVec = G4NeutronHPData::
00159          //Instance()->MakePhysicsVector((*theElementTable)[i], this);
00160          //theCrossSections->push_back(physVec);
00161       }
00162 
00163       theCrossSections->push_back(thePhysVec);
00164    }
00165 }
00166  
00167 
00168 
00169 void G4NeutronHPorLEInelasticData::DumpPhysicsTable(const G4ParticleDefinition& aP)
00170 {
00171   if(&aP!=G4Neutron::Neutron()) 
00172      throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");  
00173 //  G4cout << "G4NeutronHPorLEInelasticData::DumpPhysicsTable still to be implemented"<<G4endl;
00174 }
00175 
00176 
00177 
00178 #include "G4Nucleus.hh"
00179 #include "G4NucleiProperties.hh"
00180 #include "G4Neutron.hh"
00181 #include "G4Electron.hh"
00182 
00183 G4double G4NeutronHPorLEInelasticData::
00184 GetCrossSection(const G4DynamicParticle* aP, const G4Element*anE, G4double aT)
00185 {
00186 
00187   // G4cout << "Choice G4NeutronHPorLEInelasticData for element " << anE->GetName() << G4endl;
00188   G4double result = 0;
00189   //G4bool outOfRange;
00190   G4int index = anE->GetIndex();
00191 
00192   // prepare neutron
00193   G4double eKinetic = aP->GetKineticEnergy();
00194   G4ReactionProduct theNeutron( aP->GetDefinition() );
00195   theNeutron.SetMomentum( aP->GetMomentum() );
00196   theNeutron.SetKineticEnergy( eKinetic );
00197 
00198   // prepare thermal nucleus
00199   G4Nucleus aNuc;
00200   G4double eps = 0.0001;
00201   G4double theA = anE->GetN();
00202   G4double theZ = anE->GetZ();
00203   G4double eleMass; 
00204   eleMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theA+eps), static_cast<G4int>(theZ+eps))
00205              ) / G4Neutron::Neutron()->GetPDGMass();
00206   
00207   G4ReactionProduct boosted;
00208   G4double aXsection;
00209   
00210   // MC integration loop
00211   G4int counter = 0;
00212   G4double buffer = 0;
00213   G4int size = G4int(std::max(10., aT/60*kelvin));
00214   G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutron.GetMomentum();
00215   G4double neutronVMag = neutronVelocity.mag();
00216   while(counter == 0 || std::abs(buffer-result/std::max(1,counter)) > 0.03*buffer)
00217   {
00218     if(counter) buffer = result/counter;
00219     while (counter<size)
00220     {
00221       counter ++;
00222       G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus(eleMass, aT);
00223       boosted.Lorentz(theNeutron, aThermalNuc);
00224       G4double theEkin = boosted.GetKineticEnergy();
00225       //aXsection = (*((*theCrossSections)(index))).GetValue(theEkin, outOfRange);
00226       aXsection = theInelasticChannel[index].GetXsec( theEkin );
00227       // velocity correction.
00228       G4ThreeVector targetVelocity = 1./aThermalNuc.GetMass()*aThermalNuc.GetMomentum();
00229       aXsection *= (targetVelocity-neutronVelocity).mag()/neutronVMag;
00230       result += aXsection;
00231     }
00232     size += size;
00233   }
00234   result /= counter;
00235   //return result;
00236   return result*barn;
00237 }

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