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 #include "G4NeutronHPorLEInelastic.hh"
00038 #include "G4SystemOfUnits.hh"
00039
00040
00041 G4NeutronHPorLEInelastic::G4NeutronHPorLEInelastic()
00042 :G4HadronicInteraction("NeutronHPorLEInelastic")
00043 {
00044 SetMinEnergy(0.*eV);
00045 SetMaxEnergy(20.*MeV);
00046
00047
00048 if(!getenv("G4NEUTRONHPDATA"))
00049 throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
00050 dirName = getenv("G4NEUTRONHPDATA");
00051 G4String tString = "/Inelastic/";
00052 dirName = dirName + tString;
00053
00054 numEle = G4Element::GetNumberOfElements();
00055 theInelastic = new G4NeutronHPChannelList[numEle];
00056 unavailable_elements.clear();
00057 for (G4int i=0; i<numEle; i++)
00058 {
00059 theInelastic[i].Init( (*(G4Element::GetElementTable()))[i] , dirName );
00060 do
00061 {
00062
00063 try
00064 {
00065 theInelastic[i].Register(&theNFS, "F01");
00066 theInelastic[i].Register(&theNXFS, "F02");
00067 theInelastic[i].Register(&the2NDFS, "F03");
00068 theInelastic[i].Register(&the2NFS, "F04");
00069 theInelastic[i].Register(&the3NFS, "F05");
00070 theInelastic[i].Register(&theNAFS, "F06");
00071 theInelastic[i].Register(&theN3AFS, "F07");
00072 theInelastic[i].Register(&the2NAFS, "F08");
00073 theInelastic[i].Register(&the3NAFS, "F09");
00074 theInelastic[i].Register(&theNPFS, "F10");
00075 theInelastic[i].Register(&theN2AFS, "F11");
00076 theInelastic[i].Register(&the2N2AFS, "F12");
00077 theInelastic[i].Register(&theNDFS, "F13");
00078 theInelastic[i].Register(&theNTFS, "F14");
00079 theInelastic[i].Register(&theNHe3FS, "F15");
00080 theInelastic[i].Register(&theND2AFS, "F16");
00081 theInelastic[i].Register(&theNT2AFS, "F17");
00082 theInelastic[i].Register(&the4NFS, "F18");
00083 theInelastic[i].Register(&the2NPFS, "F19");
00084 theInelastic[i].Register(&the3NPFS, "F20");
00085 theInelastic[i].Register(&theN2PFS, "F21");
00086 theInelastic[i].Register(&theNPAFS, "F22");
00087 theInelastic[i].Register(&thePFS, "F23");
00088 theInelastic[i].Register(&theDFS, "F24");
00089 theInelastic[i].Register(&theTFS, "F25");
00090 theInelastic[i].Register(&theHe3FS, "F26");
00091 theInelastic[i].Register(&theAFS, "F27");
00092 theInelastic[i].Register(&the2AFS, "F28");
00093 theInelastic[i].Register(&the3AFS, "F29");
00094 theInelastic[i].Register(&the2PFS, "F30");
00095 theInelastic[i].Register(&thePAFS, "F31");
00096 theInelastic[i].Register(&theD2AFS, "F32");
00097 theInelastic[i].Register(&theT2AFS, "F33");
00098 theInelastic[i].Register(&thePDFS, "F34");
00099 theInelastic[i].Register(&thePTFS, "F35");
00100 theInelastic[i].Register(&theDAFS, "F36");
00101 }
00102
00103 catch ( G4HadronicException )
00104 {
00105 unavailable_elements.insert ( (*(G4Element::GetElementTable()))[i]->GetName() );
00106 }
00107 theInelastic[i].RestartRegistration();
00108 }
00109 while( !theInelastic[i].HasDataInAnyFinalState());
00110
00111 }
00112
00113
00114 if ( unavailable_elements.size() > 0 )
00115 {
00116 std::set< G4String>::iterator it;
00117 G4cout << "HP Inelastic data are not available for thess elements "<< G4endl;
00118 for ( it = unavailable_elements.begin() ; it != unavailable_elements.end() ; it++ )
00119 G4cout << *it << G4endl;
00120 G4cout << "Low Energy Parameterization Models will be used."<< G4endl;
00121 }
00122
00123 createXSectionDataSet();
00124 }
00125
00126 G4NeutronHPorLEInelastic::~G4NeutronHPorLEInelastic()
00127 {
00128 delete [] theInelastic;
00129 delete theDataSet;
00130 }
00131
00132 #include "G4NeutronHPThermalBoost.hh"
00133
00134 G4HadFinalState * G4NeutronHPorLEInelastic::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& )
00135 {
00136 G4int it=0;
00137 const G4Material * theMaterial = aTrack.GetMaterial();
00138 G4int n = theMaterial->GetNumberOfElements();
00139 G4int index = theMaterial->GetElement(0)->GetIndex();
00140 if(n!=1)
00141 {
00142 G4int i;
00143 xSec = new G4double[n];
00144 G4double sum=0;
00145 const G4double * NumAtomsPerVolume = theMaterial->GetVecNbOfAtomsPerVolume();
00146 G4double rWeight;
00147 G4NeutronHPThermalBoost aThermalE;
00148 for (i=0; i<n; i++)
00149 {
00150 index = theMaterial->GetElement(i)->GetIndex();
00151 rWeight = NumAtomsPerVolume[i];
00152 G4double x = aThermalE.GetThermalEnergy(aTrack, theMaterial->GetElement(i), theMaterial->GetTemperature());
00153
00154
00155
00156
00157 xSec[i] = theInelastic[index].GetXsec(x);
00158
00159 xSec[i] *= rWeight;
00160 sum+=xSec[i];
00161 }
00162 G4double random = G4UniformRand();
00163 G4double running = 0;
00164 for (i=0; i<n; i++)
00165 {
00166 running += xSec[i];
00167 index = theMaterial->GetElement(i)->GetIndex();
00168 it = i;
00169 if(random<=running/sum) break;
00170 }
00171 delete [] xSec;
00172
00173 }
00174
00175 return theInelastic[index].ApplyYourself( theMaterial->GetElement(it) , aTrack );
00176 }
00177
00178
00179
00180 G4bool G4NeutronHPorLEInelastic::IsThisElementOK( G4String name )
00181 {
00182 if ( unavailable_elements.find( name ) == unavailable_elements.end() )
00183 return TRUE;
00184 else
00185 return FALSE;
00186 }
00187
00188
00189
00190 void G4NeutronHPorLEInelastic::createXSectionDataSet()
00191 {
00192 theDataSet = new G4NeutronHPorLEInelasticData ( theInelastic , &unavailable_elements );
00193 }
00194 const std::pair<G4double, G4double> G4NeutronHPorLEInelastic::GetFatalEnergyCheckLevels() const
00195 {
00196
00197 return std::pair<G4double, G4double>(10*perCent,DBL_MAX);
00198 }