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 #include "G4NeutronHPorLElastic.hh"
00039 #include "G4SystemOfUnits.hh"
00040 #include "G4NeutronHPElasticFS.hh"
00041
00042 G4NeutronHPorLElastic::G4NeutronHPorLElastic()
00043 :G4HadronicInteraction("NeutronHPorLElastic")
00044 {
00045 overrideSuspension = false;
00046 G4NeutronHPElasticFS * theFS = new G4NeutronHPElasticFS;
00047 if(!getenv("G4NEUTRONHPDATA"))
00048 throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
00049 dirName = getenv("G4NEUTRONHPDATA");
00050 G4String tString = "/Elastic/";
00051 dirName = dirName + tString;
00052
00053 numEle = G4Element::GetNumberOfElements();
00054 theElastic = new G4NeutronHPChannel[numEle];
00055 unavailable_elements.clear();
00056 for (G4int i=0; i<numEle; i++)
00057 {
00058 theElastic[i].Init((*(G4Element::GetElementTable()))[i], dirName);
00059
00060 try { while(!theElastic[i].Register(theFS)) ; }
00061 catch ( G4HadronicException )
00062 {
00063 unavailable_elements.insert ( (*(G4Element::GetElementTable()))[i]->GetName() );
00064 }
00065 }
00066 delete theFS;
00067 SetMinEnergy(0.*eV);
00068 SetMaxEnergy(20.*MeV);
00069 if ( unavailable_elements.size() > 0 )
00070 {
00071 std::set< G4String>::iterator it;
00072 G4cout << "HP Elastic data are not available for thess elements "<< G4endl;
00073 for ( it = unavailable_elements.begin() ; it != unavailable_elements.end() ; it++ )
00074 G4cout << *it << G4endl;
00075 G4cout << "Low Energy Parameterization Models will be used."<< G4endl;
00076 }
00077
00078 createXSectionDataSet();
00079 }
00080
00081 G4NeutronHPorLElastic::~G4NeutronHPorLElastic()
00082 {
00083 delete [] theElastic;
00084 delete theDataSet;
00085 }
00086
00087 #include "G4NeutronHPThermalBoost.hh"
00088
00089 G4HadFinalState * G4NeutronHPorLElastic::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& )
00090 {
00091 const G4Material * theMaterial = aTrack.GetMaterial();
00092 G4int n = theMaterial->GetNumberOfElements();
00093 G4int index = theMaterial->GetElement(0)->GetIndex();
00094 if(n!=1)
00095 {
00096 G4int i;
00097 xSec = new G4double[n];
00098 G4double sum=0;
00099 const G4double * NumAtomsPerVolume = theMaterial->GetVecNbOfAtomsPerVolume();
00100 G4double rWeight;
00101 G4NeutronHPThermalBoost aThermalE;
00102 for (i=0; i<n; i++)
00103 {
00104 index = theMaterial->GetElement(i)->GetIndex();
00105 rWeight = NumAtomsPerVolume[i];
00106 G4double x = aThermalE.GetThermalEnergy(aTrack, theMaterial->GetElement(i), theMaterial->GetTemperature());
00107
00108
00109
00110
00111 xSec[i] = theElastic[index].GetXsec(x);
00112
00113 xSec[i] *= rWeight;
00114 sum+=xSec[i];
00115 }
00116 G4double random = G4UniformRand();
00117 G4double running = 0;
00118 for (i=0; i<n; i++)
00119 {
00120 running += xSec[i];
00121 index = theMaterial->GetElement(i)->GetIndex();
00122 if(random<=running/sum) break;
00123 }
00124 delete [] xSec;
00125
00126 }
00127 G4HadFinalState* finalState = theElastic[index].ApplyYourself(aTrack);
00128 if (overrideSuspension) finalState->SetStatusChange(isAlive);
00129 return finalState;
00130 }
00131
00132
00133
00134 G4bool G4NeutronHPorLElastic::IsThisElementOK( G4String name )
00135 {
00136 if ( unavailable_elements.find( name ) == unavailable_elements.end() )
00137 return TRUE;
00138 else
00139 return FALSE;
00140 }
00141
00142
00143
00144 void G4NeutronHPorLElastic::createXSectionDataSet()
00145 {
00146 theDataSet = new G4NeutronHPorLElasticData ( theElastic , &unavailable_elements );
00147 }
00148
00149 const std::pair<G4double, G4double> G4NeutronHPorLElastic::GetFatalEnergyCheckLevels() const
00150 {
00151
00152 return std::pair<G4double, G4double>(10*perCent,DBL_MAX);
00153 }