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
00040 #include "G4NeutronHPInelastic.hh"
00041 #include "G4SystemOfUnits.hh"
00042
00043 #include "G4NeutronHPManager.hh"
00044
00045 G4NeutronHPInelastic::G4NeutronHPInelastic()
00046 :G4HadronicInteraction("NeutronHPInelastic")
00047 {
00048 SetMinEnergy( 0.0 );
00049 SetMaxEnergy( 20.*MeV );
00050
00051 G4int istatus = system("echo $G4NEUTRONHPDATA");
00052 if ( istatus < 0 )
00053 {
00054 G4cout << "Warning! system(\"echo $G4NEUTRONHPDATA\") returns error value at G4NeutronHPInelastic" << G4endl;
00055 }
00056
00057
00058 if(!getenv("G4NEUTRONHPDATA"))
00059 throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
00060 dirName = getenv("G4NEUTRONHPDATA");
00061 G4String tString = "/Inelastic";
00062 dirName = dirName + tString;
00063 numEle = G4Element::GetNumberOfElements();
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128 for (G4int i=0; i<numEle; i++)
00129 {
00130 theInelastic.push_back( new G4NeutronHPChannelList );
00131 (*theInelastic[i]).Init((*(G4Element::GetElementTable()))[i], dirName);
00132 G4int itry = 0;
00133 do
00134 {
00135 (*theInelastic[i]).Register(&theNFS, "F01");
00136 (*theInelastic[i]).Register(&theNXFS, "F02");
00137 (*theInelastic[i]).Register(&the2NDFS, "F03");
00138 (*theInelastic[i]).Register(&the2NFS, "F04");
00139 (*theInelastic[i]).Register(&the3NFS, "F05");
00140 (*theInelastic[i]).Register(&theNAFS, "F06");
00141 (*theInelastic[i]).Register(&theN3AFS, "F07");
00142 (*theInelastic[i]).Register(&the2NAFS, "F08");
00143 (*theInelastic[i]).Register(&the3NAFS, "F09");
00144 (*theInelastic[i]).Register(&theNPFS, "F10");
00145 (*theInelastic[i]).Register(&theN2AFS, "F11");
00146 (*theInelastic[i]).Register(&the2N2AFS, "F12");
00147 (*theInelastic[i]).Register(&theNDFS, "F13");
00148 (*theInelastic[i]).Register(&theNTFS, "F14");
00149 (*theInelastic[i]).Register(&theNHe3FS, "F15");
00150 (*theInelastic[i]).Register(&theND2AFS, "F16");
00151 (*theInelastic[i]).Register(&theNT2AFS, "F17");
00152 (*theInelastic[i]).Register(&the4NFS, "F18");
00153 (*theInelastic[i]).Register(&the2NPFS, "F19");
00154 (*theInelastic[i]).Register(&the3NPFS, "F20");
00155 (*theInelastic[i]).Register(&theN2PFS, "F21");
00156 (*theInelastic[i]).Register(&theNPAFS, "F22");
00157 (*theInelastic[i]).Register(&thePFS, "F23");
00158 (*theInelastic[i]).Register(&theDFS, "F24");
00159 (*theInelastic[i]).Register(&theTFS, "F25");
00160 (*theInelastic[i]).Register(&theHe3FS, "F26");
00161 (*theInelastic[i]).Register(&theAFS, "F27");
00162 (*theInelastic[i]).Register(&the2AFS, "F28");
00163 (*theInelastic[i]).Register(&the3AFS, "F29");
00164 (*theInelastic[i]).Register(&the2PFS, "F30");
00165 (*theInelastic[i]).Register(&thePAFS, "F31");
00166 (*theInelastic[i]).Register(&theD2AFS, "F32");
00167 (*theInelastic[i]).Register(&theT2AFS, "F33");
00168 (*theInelastic[i]).Register(&thePDFS, "F34");
00169 (*theInelastic[i]).Register(&thePTFS, "F35");
00170 (*theInelastic[i]).Register(&theDAFS, "F36");
00171 (*theInelastic[i]).RestartRegistration();
00172 itry++;
00173 }
00174 while( !(*theInelastic[i]).HasDataInAnyFinalState() && itry < 6 );
00175
00176
00177 if ( itry == 6 )
00178 {
00179
00180 G4bool exceptional = false;
00181 if ( (*(G4Element::GetElementTable()))[i]->GetNumberOfIsotopes() == 1 )
00182 {
00183 if ( (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetZ() == 1 && (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetN() == 1 ) exceptional = true;
00184 }
00185 if ( !exceptional ) throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this element");
00186 }
00187
00188 }
00189 }
00190
00191 G4NeutronHPInelastic::~G4NeutronHPInelastic()
00192 {
00193
00194 for ( std::vector<G4NeutronHPChannelList*>::iterator
00195 it = theInelastic.begin() ; it != theInelastic.end() ; it++ )
00196 {
00197 delete *it;
00198 }
00199 theInelastic.clear();
00200 }
00201
00202 #include "G4NeutronHPThermalBoost.hh"
00203
00204 G4HadFinalState * G4NeutronHPInelastic::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aNucleus )
00205 {
00206 if ( numEle < (G4int)G4Element::GetNumberOfElements() ) addChannelForNewElement();
00207 G4NeutronHPManager::GetInstance()->OpenReactionWhiteBoard();
00208 const G4Material * theMaterial = aTrack.GetMaterial();
00209 G4int n = theMaterial->GetNumberOfElements();
00210 G4int index = theMaterial->GetElement(0)->GetIndex();
00211 G4int it=0;
00212 if(n!=1)
00213 {
00214 xSec = new G4double[n];
00215 G4double sum=0;
00216 G4int i;
00217 const G4double * NumAtomsPerVolume = theMaterial->GetVecNbOfAtomsPerVolume();
00218 G4double rWeight;
00219 G4NeutronHPThermalBoost aThermalE;
00220 for (i=0; i<n; i++)
00221 {
00222 index = theMaterial->GetElement(i)->GetIndex();
00223 rWeight = NumAtomsPerVolume[i];
00224
00225 xSec[i] = (*theInelastic[index]).GetXsec(aThermalE.GetThermalEnergy(aTrack,
00226 theMaterial->GetElement(i),
00227 theMaterial->GetTemperature()));
00228 xSec[i] *= rWeight;
00229 sum+=xSec[i];
00230 }
00231 G4double random = G4UniformRand();
00232 G4double running = 0;
00233 for (i=0; i<n; i++)
00234 {
00235 running += xSec[i];
00236 index = theMaterial->GetElement(i)->GetIndex();
00237 it = i;
00238
00239 if( sum == 0 || random<=running/sum) break;
00240 }
00241 delete [] xSec;
00242 }
00243
00244
00245 G4HadFinalState* result = (*theInelastic[index]).ApplyYourself(theMaterial->GetElement(it), aTrack);
00246
00247 aNucleus.SetParameters(G4NeutronHPManager::GetInstance()->GetReactionWhiteBoard()->GetTargA(),G4NeutronHPManager::GetInstance()->GetReactionWhiteBoard()->GetTargZ());
00248 G4NeutronHPManager::GetInstance()->CloseReactionWhiteBoard();
00249 return result;
00250 }
00251
00252 const std::pair<G4double, G4double> G4NeutronHPInelastic::GetFatalEnergyCheckLevels() const
00253 {
00254
00255
00256
00257
00258 return std::pair<G4double, G4double>(10*perCent,DBL_MAX);
00259 }
00260
00261 void G4NeutronHPInelastic::addChannelForNewElement()
00262 {
00263 for ( G4int i = numEle ; i < (G4int)G4Element::GetNumberOfElements() ; i++ )
00264 {
00265 G4cout << "G4NeutronHPInelastic Prepairing Data for the new element of " << (*(G4Element::GetElementTable()))[i]->GetName() << G4endl;
00266
00267 theInelastic.push_back( new G4NeutronHPChannelList );
00268 (*theInelastic[i]).Init((*(G4Element::GetElementTable()))[i], dirName);
00269 G4int itry = 0;
00270 do
00271 {
00272 (*theInelastic[i]).Register(&theNFS, "F01");
00273 (*theInelastic[i]).Register(&theNXFS, "F02");
00274 (*theInelastic[i]).Register(&the2NDFS, "F03");
00275 (*theInelastic[i]).Register(&the2NFS, "F04");
00276 (*theInelastic[i]).Register(&the3NFS, "F05");
00277 (*theInelastic[i]).Register(&theNAFS, "F06");
00278 (*theInelastic[i]).Register(&theN3AFS, "F07");
00279 (*theInelastic[i]).Register(&the2NAFS, "F08");
00280 (*theInelastic[i]).Register(&the3NAFS, "F09");
00281 (*theInelastic[i]).Register(&theNPFS, "F10");
00282 (*theInelastic[i]).Register(&theN2AFS, "F11");
00283 (*theInelastic[i]).Register(&the2N2AFS, "F12");
00284 (*theInelastic[i]).Register(&theNDFS, "F13");
00285 (*theInelastic[i]).Register(&theNTFS, "F14");
00286 (*theInelastic[i]).Register(&theNHe3FS, "F15");
00287 (*theInelastic[i]).Register(&theND2AFS, "F16");
00288 (*theInelastic[i]).Register(&theNT2AFS, "F17");
00289 (*theInelastic[i]).Register(&the4NFS, "F18");
00290 (*theInelastic[i]).Register(&the2NPFS, "F19");
00291 (*theInelastic[i]).Register(&the3NPFS, "F20");
00292 (*theInelastic[i]).Register(&theN2PFS, "F21");
00293 (*theInelastic[i]).Register(&theNPAFS, "F22");
00294 (*theInelastic[i]).Register(&thePFS, "F23");
00295 (*theInelastic[i]).Register(&theDFS, "F24");
00296 (*theInelastic[i]).Register(&theTFS, "F25");
00297 (*theInelastic[i]).Register(&theHe3FS, "F26");
00298 (*theInelastic[i]).Register(&theAFS, "F27");
00299 (*theInelastic[i]).Register(&the2AFS, "F28");
00300 (*theInelastic[i]).Register(&the3AFS, "F29");
00301 (*theInelastic[i]).Register(&the2PFS, "F30");
00302 (*theInelastic[i]).Register(&thePAFS, "F31");
00303 (*theInelastic[i]).Register(&theD2AFS, "F32");
00304 (*theInelastic[i]).Register(&theT2AFS, "F33");
00305 (*theInelastic[i]).Register(&thePDFS, "F34");
00306 (*theInelastic[i]).Register(&thePTFS, "F35");
00307 (*theInelastic[i]).Register(&theDAFS, "F36");
00308 (*theInelastic[i]).RestartRegistration();
00309 itry++;
00310 }
00311 while( !(*theInelastic[i]).HasDataInAnyFinalState() && itry < 6 );
00312
00313
00314 if ( itry == 6 )
00315 {
00316
00317 G4bool exceptional = false;
00318 if ( (*(G4Element::GetElementTable()))[i]->GetNumberOfIsotopes() == 1 )
00319 {
00320 if ( (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetZ() == 1 && (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetN() == 1 ) exceptional = true;
00321 }
00322 if ( !exceptional ) throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this element");
00323 }
00324 }
00325
00326 numEle = (G4int)G4Element::GetNumberOfElements();
00327 }