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 #include "G4NeutronHPElastic.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "G4NeutronHPElasticFS.hh"
00036 #include "G4NeutronHPManager.hh"
00037
00038 G4NeutronHPElastic::G4NeutronHPElastic()
00039 :G4HadronicInteraction("NeutronHPElastic")
00040 {
00041 overrideSuspension = false;
00042 G4NeutronHPElasticFS * theFS = new G4NeutronHPElasticFS;
00043 if(!getenv("G4NEUTRONHPDATA"))
00044 throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
00045 dirName = getenv("G4NEUTRONHPDATA");
00046 G4String tString = "/Elastic";
00047 dirName = dirName + tString;
00048
00049 numEle = G4Element::GetNumberOfElements();
00050
00051
00052
00053
00054
00055
00056 for ( G4int i = 0 ; i < numEle ; i++ )
00057 {
00058 theElastic.push_back( new G4NeutronHPChannel );
00059 (*theElastic[i]).Init((*(G4Element::GetElementTable()))[i], dirName);
00060 while(!(*theElastic[i]).Register(theFS)) ;
00061 }
00062 delete theFS;
00063 SetMinEnergy(0.*eV);
00064 SetMaxEnergy(20.*MeV);
00065 }
00066
00067 G4NeutronHPElastic::~G4NeutronHPElastic()
00068 {
00069
00070 for ( std::vector<G4NeutronHPChannel*>::iterator
00071 it = theElastic.begin() ; it != theElastic.end() ; it++ )
00072 {
00073 delete *it;
00074 }
00075 theElastic.clear();
00076 }
00077
00078 #include "G4NeutronHPThermalBoost.hh"
00079
00080 G4HadFinalState * G4NeutronHPElastic::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aNucleus )
00081 {
00082
00083 if ( numEle < (G4int)G4Element::GetNumberOfElements() ) addChannelForNewElement();
00084
00085 G4NeutronHPManager::GetInstance()->OpenReactionWhiteBoard();
00086 const G4Material * theMaterial = aTrack.GetMaterial();
00087 G4int n = theMaterial->GetNumberOfElements();
00088 G4int index = theMaterial->GetElement(0)->GetIndex();
00089 if(n!=1)
00090 {
00091 G4int i;
00092 xSec = new G4double[n];
00093 G4double sum=0;
00094 const G4double * NumAtomsPerVolume = theMaterial->GetVecNbOfAtomsPerVolume();
00095 G4double rWeight;
00096 G4NeutronHPThermalBoost aThermalE;
00097 for (i=0; i<n; i++)
00098 {
00099 index = theMaterial->GetElement(i)->GetIndex();
00100 rWeight = NumAtomsPerVolume[i];
00101
00102 xSec[i] = (*theElastic[index]).GetXsec(aThermalE.GetThermalEnergy(aTrack,
00103 theMaterial->GetElement(i),
00104 theMaterial->GetTemperature()));
00105 xSec[i] *= rWeight;
00106 sum+=xSec[i];
00107 }
00108 G4double random = G4UniformRand();
00109 G4double running = 0;
00110 for (i=0; i<n; i++)
00111 {
00112 running += xSec[i];
00113 index = theMaterial->GetElement(i)->GetIndex();
00114
00115 if( sum == 0 || random <= running/sum ) break;
00116 }
00117 delete [] xSec;
00118
00119 }
00120
00121 G4HadFinalState* finalState = (*theElastic[index]).ApplyYourself(aTrack);
00122 if (overrideSuspension) finalState->SetStatusChange(isAlive);
00123 aNucleus.SetParameters(G4NeutronHPManager::GetInstance()->GetReactionWhiteBoard()->GetTargA(),G4NeutronHPManager::GetInstance()->GetReactionWhiteBoard()->GetTargZ());
00124 G4NeutronHPManager::GetInstance()->CloseReactionWhiteBoard();
00125 return finalState;
00126 }
00127
00128 const std::pair<G4double, G4double> G4NeutronHPElastic::GetFatalEnergyCheckLevels() const
00129 {
00130
00131 return std::pair<G4double, G4double>(10*perCent,DBL_MAX);
00132 }
00133
00134 void G4NeutronHPElastic::addChannelForNewElement()
00135 {
00136 G4NeutronHPElasticFS* theFS = new G4NeutronHPElasticFS;
00137 for ( G4int i = numEle ; i < (G4int)G4Element::GetNumberOfElements() ; i++ )
00138 {
00139 G4cout << "G4NeutronHPElastic Prepairing Data for the new element of " << (*(G4Element::GetElementTable()))[i]->GetName() << G4endl;
00140 theElastic.push_back( new G4NeutronHPChannel );
00141 (*theElastic[i]).Init((*(G4Element::GetElementTable()))[i], dirName);
00142 while(!(*theElastic[i]).Register(theFS)) ;
00143 }
00144 delete theFS;
00145 numEle = (G4int)G4Element::GetNumberOfElements();
00146 }