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