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
00041
00042
00043
00044 #include "G4NeutronHPThermalScattering.hh"
00045 #include "G4SystemOfUnits.hh"
00046 #include "G4Neutron.hh"
00047 #include "G4ElementTable.hh"
00048
00049 G4NeutronHPThermalScattering::G4NeutronHPThermalScattering()
00050 :G4HadronicInteraction("NeutronHPThermalScattering")
00051 {
00052
00053 theHPElastic = new G4NeutronHPElastic();
00054
00055 SetMinEnergy( 0.*eV );
00056 SetMaxEnergy( 4*eV );
00057 theXSection = new G4NeutronHPThermalScatteringData();
00058 theXSection->BuildPhysicsTable( *(G4Neutron::Neutron()) );
00059
00060 buildPhysicsTable();
00061
00062 }
00063
00064
00065
00066 G4NeutronHPThermalScattering::~G4NeutronHPThermalScattering()
00067 {
00068
00069 for ( std::map < G4int , std::map < G4double , std::vector < E_isoAng* >* >* >::iterator it = incoherentFSs.begin() ; it != incoherentFSs.end() ; it++ )
00070 {
00071 std::map < G4double , std::vector < E_isoAng* >* >::iterator itt;
00072 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
00073 {
00074 std::vector< E_isoAng* >::iterator ittt;
00075 for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
00076 {
00077 delete *ittt;
00078 }
00079 delete itt->second;
00080 }
00081 delete it->second;
00082 }
00083
00084 for ( std::map < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* >::iterator it = coherentFSs.begin() ; it != coherentFSs.end() ; it++ )
00085 {
00086 std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >::iterator itt;
00087 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
00088 {
00089 std::vector < std::pair< G4double , G4double >* >::iterator ittt;
00090 for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
00091 {
00092 delete *ittt;
00093 }
00094 delete itt->second;
00095 }
00096 delete it->second;
00097 }
00098
00099 for ( std::map < G4int , std::map < G4double , std::vector < E_P_E_isoAng* >* >* >::iterator it = inelasticFSs.begin() ; it != inelasticFSs.end() ; it++ )
00100 {
00101 std::map < G4double , std::vector < E_P_E_isoAng* >* >::iterator itt;
00102 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
00103 {
00104 std::vector < E_P_E_isoAng* >::iterator ittt;
00105 for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
00106 {
00107 std::vector < E_isoAng* >::iterator it4;
00108 for ( it4 = (*ittt)->vE_isoAngle.begin() ; it4 != (*ittt)->vE_isoAngle.end() ; it4++ )
00109 {
00110 delete *it4;
00111 }
00112 delete *ittt;
00113 }
00114 delete itt->second;
00115 }
00116 delete it->second;
00117 }
00118
00119 delete theHPElastic;
00120 delete theXSection;
00121 }
00122
00123
00124
00125 std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* G4NeutronHPThermalScattering::readACoherentFSDATA( G4String name )
00126 {
00127
00128 std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* aCoherentFSDATA = new std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >;
00129
00130 std::ifstream theChannel( name.c_str() );
00131
00132 std::vector< G4double > vBraggE;
00133
00134 G4int dummy;
00135 while ( theChannel >> dummy )
00136 {
00137 theChannel >> dummy;
00138 G4double temp;
00139 theChannel >> temp;
00140 std::vector < std::pair< G4double , G4double >* >* anBragE_P = new std::vector < std::pair< G4double , G4double >* >;
00141
00142 G4int n;
00143 theChannel >> n;
00144 for ( G4int i = 0 ; i < n ; i++ )
00145 {
00146 G4double Ei;
00147 G4double Pi;
00148 if ( aCoherentFSDATA->size() == 0 )
00149 {
00150 theChannel >> Ei;
00151 vBraggE.push_back( Ei );
00152 }
00153 else
00154 {
00155 Ei = vBraggE[ i ];
00156 }
00157 theChannel >> Pi;
00158 anBragE_P->push_back ( new std::pair < G4double , G4double > ( Ei , Pi ) );
00159
00160 }
00161 aCoherentFSDATA->insert ( std::pair < G4double , std::vector < std::pair< G4double , G4double >* >* > ( temp , anBragE_P ) );
00162 }
00163
00164 return aCoherentFSDATA;
00165 }
00166
00167
00168
00169 std::map < G4double , std::vector < E_P_E_isoAng* >* >* G4NeutronHPThermalScattering::readAnInelasticFSDATA ( G4String name )
00170 {
00171 std::map < G4double , std::vector < E_P_E_isoAng* >* >* anT_E_P_E_isoAng = new std::map < G4double , std::vector < E_P_E_isoAng* >* >;
00172
00173 std::ifstream theChannel( name.c_str() );
00174
00175 G4int dummy;
00176 while ( theChannel >> dummy )
00177 {
00178 theChannel >> dummy;
00179 G4double temp;
00180 theChannel >> temp;
00181 std::vector < E_P_E_isoAng* >* vE_P_E_isoAng = new std::vector < E_P_E_isoAng* >;
00182 G4int n;
00183 theChannel >> n;
00184 for ( G4int i = 0 ; i < n ; i++ )
00185 {
00186 vE_P_E_isoAng->push_back ( readAnE_P_E_isoAng ( &theChannel ) );
00187 }
00188 anT_E_P_E_isoAng->insert ( std::pair < G4double , std::vector < E_P_E_isoAng* >* > ( temp , vE_P_E_isoAng ) );
00189 }
00190 theChannel.close();
00191
00192 return anT_E_P_E_isoAng;
00193 }
00194
00195
00196
00197 E_P_E_isoAng* G4NeutronHPThermalScattering::readAnE_P_E_isoAng( std::ifstream* file )
00198 {
00199 E_P_E_isoAng* aData = new E_P_E_isoAng;
00200
00201 G4double dummy;
00202 G4double energy;
00203 G4int nep , nl;
00204 *file >> dummy;
00205 *file >> energy;
00206 aData->energy = energy*eV;
00207 *file >> dummy;
00208 *file >> dummy;
00209 *file >> nep;
00210 *file >> nl;
00211 aData->n = nep/nl;
00212 for ( G4int i = 0 ; i < aData->n ; i++ )
00213 {
00214 G4double prob;
00215 E_isoAng* anE_isoAng = new E_isoAng;
00216 aData->vE_isoAngle.push_back( anE_isoAng );
00217 *file >> energy;
00218 anE_isoAng->energy = energy*eV;
00219 anE_isoAng->n = nl - 2;
00220 anE_isoAng->isoAngle.resize( anE_isoAng->n );
00221 *file >> prob;
00222 aData->prob.push_back( prob );
00223
00224 for ( G4int j = 0 ; j < anE_isoAng->n ; j++ )
00225 {
00226 G4double x;
00227 *file >> x;
00228 anE_isoAng->isoAngle[j] = x ;
00229
00230 }
00231 }
00232
00233
00234 G4double total = 0;
00235 for ( G4int i = 0 ; i < aData->n - 1 ; i++ )
00236 {
00237 G4double E_L = aData->vE_isoAngle[i]->energy/eV;
00238 G4double E_H = aData->vE_isoAngle[i+1]->energy/eV;
00239 G4double dE = E_H - E_L;
00240 total += ( ( aData->prob[i] ) * dE );
00241 }
00242 aData->sum_of_probXdEs = total;
00243
00244 return aData;
00245 }
00246
00247
00248
00249 std::map < G4double , std::vector < E_isoAng* >* >* G4NeutronHPThermalScattering::readAnIncoherentFSDATA ( G4String name )
00250 {
00251 std::map < G4double , std::vector < E_isoAng* >* >* T_E = new std::map < G4double , std::vector < E_isoAng* >* >;
00252
00253 std::ifstream theChannel( name.c_str() );
00254
00255 G4int dummy;
00256 while ( theChannel >> dummy )
00257 {
00258 theChannel >> dummy;
00259 G4double temp;
00260 theChannel >> temp;
00261 std::vector < E_isoAng* >* vE_isoAng = new std::vector < E_isoAng* >;
00262 G4int n;
00263 theChannel >> n;
00264 for ( G4int i = 0 ; i < n ; i++ )
00265 vE_isoAng->push_back ( readAnE_isoAng( &theChannel ) );
00266 T_E->insert ( std::pair < G4double , std::vector < E_isoAng* >* > ( temp , vE_isoAng ) );
00267 }
00268 theChannel.close();
00269
00270 return T_E;
00271 }
00272
00273
00274
00275 E_isoAng* G4NeutronHPThermalScattering::readAnE_isoAng( std::ifstream* file )
00276 {
00277 E_isoAng* aData = new E_isoAng;
00278
00279 G4double dummy;
00280 G4double energy;
00281 G4int n;
00282 *file >> dummy;
00283 *file >> energy;
00284 *file >> dummy;
00285 *file >> dummy;
00286 *file >> n;
00287 *file >> dummy;
00288 aData->energy = energy*eV;
00289 aData->n = n-2;
00290 aData->isoAngle.resize( n );
00291
00292 *file >> dummy;
00293 *file >> dummy;
00294 for ( G4int i = 0 ; i < aData->n ; i++ )
00295 *file >> aData->isoAngle[i];
00296
00297 return aData;
00298 }
00299
00300
00301
00302 G4HadFinalState* G4NeutronHPThermalScattering::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aNucleus )
00303 {
00304
00305
00306
00307 const G4Material * theMaterial = aTrack.GetMaterial();
00308 G4double aTemp = theMaterial->GetTemperature();
00309 G4int n = theMaterial->GetNumberOfElements();
00310
00311
00312 G4bool findThermalElement = false;
00313 G4int ielement;
00314 const G4Element* theElement = NULL;
00315 for ( G4int i = 0; i < n ; i++ )
00316 {
00317 theElement = theMaterial->GetElement(i);
00318
00319 if ( aNucleus.GetZ_asInt() == (G4int)(theElement->GetZ() + 0.5 ) )
00320 {
00321
00322 if ( getTS_ID( NULL , theElement ) != -1 )
00323 {
00324 ielement = getTS_ID( NULL , theElement );
00325 findThermalElement = true;
00326 break;
00327 }
00328 else if ( getTS_ID( theMaterial , theElement ) != -1 )
00329 {
00330 ielement = getTS_ID( theMaterial , theElement );
00331 findThermalElement = true;
00332 break;
00333 }
00334 }
00335 }
00336
00337 if ( findThermalElement == true )
00338 {
00339
00340
00341
00342 G4ParticleDefinition* pd = const_cast< G4ParticleDefinition* >( aTrack.GetDefinition() );
00343 G4DynamicParticle* dp = new G4DynamicParticle ( pd , aTrack.Get4Momentum() );
00344 G4double total = theXSection->GetCrossSection( dp , theElement , theMaterial );
00345 G4double inelastic = theXSection->GetInelasticCrossSection( dp , theElement , theMaterial );
00346
00347
00348 G4double random = G4UniformRand();
00349 if ( random <= inelastic/total )
00350 {
00351
00352
00353
00354 std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator it;
00355 std::vector<G4double> v_temp;
00356 v_temp.clear();
00357 for ( it = inelasticFSs.find( ielement )->second->begin() ; it != inelasticFSs.find( ielement )->second->end() ; it++ )
00358 {
00359 v_temp.push_back( it->first );
00360 }
00361
00362
00363 std::pair < G4double , G4double > tempLH = find_LH ( aTemp , &v_temp );
00364
00365
00366
00367 std::vector< E_P_E_isoAng* >* vNEP_EPM_TL = 0;
00368 std::vector< E_P_E_isoAng* >* vNEP_EPM_TH = 0;
00369
00370 if ( tempLH.first != 0.0 && tempLH.second != 0.0 )
00371 {
00372 vNEP_EPM_TL = inelasticFSs.find( ielement )->second->find ( tempLH.first/kelvin )->second;
00373 vNEP_EPM_TH = inelasticFSs.find( ielement )->second->find ( tempLH.second/kelvin )->second;
00374 }
00375 else if ( tempLH.first == 0.0 )
00376 {
00377 std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator itm;
00378 itm = inelasticFSs.find( ielement )->second->begin();
00379 vNEP_EPM_TL = itm->second;
00380 itm++;
00381 vNEP_EPM_TH = itm->second;
00382 }
00383 else if ( tempLH.second == 0.0 )
00384 {
00385 std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator itm;
00386 itm = inelasticFSs.find( ielement )->second->end();
00387 itm--;
00388 vNEP_EPM_TH = itm->second;
00389 itm--;
00390 vNEP_EPM_TL = itm->second;
00391 }
00392
00393
00394
00395 G4double rand_for_sE = G4UniformRand();
00396
00397 std::pair< G4double , E_isoAng > TL = create_sE_and_EPM_from_pE_and_vE_P_E_isoAng ( rand_for_sE , aTrack.GetKineticEnergy() , vNEP_EPM_TL );
00398 std::pair< G4double , E_isoAng > TH = create_sE_and_EPM_from_pE_and_vE_P_E_isoAng ( rand_for_sE , aTrack.GetKineticEnergy() , vNEP_EPM_TH );
00399
00400 G4double sE;
00401 sE = get_linear_interpolated ( aTemp , std::pair < G4double , G4double > ( tempLH.first , TL.first ) , std::pair < G4double , G4double > ( tempLH.second , TH.first ) );
00402 E_isoAng anE_isoAng;
00403 if ( TL.second.n == TH.second.n )
00404 {
00405 anE_isoAng.energy = sE;
00406 anE_isoAng.n = TL.second.n;
00407 for ( G4int i=0 ; i < anE_isoAng.n ; i++ )
00408 {
00409 G4double angle;
00410 angle = get_linear_interpolated ( aTemp , std::pair< G4double , G4double > ( tempLH.first , TL.second.isoAngle[ i ] ) , std::pair< G4double , G4double > ( tempLH.second , TH.second.isoAngle[ i ] ) );
00411 anE_isoAng.isoAngle.push_back( angle );
00412 }
00413 }
00414 else
00415 {
00416 std::cout << "Do not Suuport yet." << std::endl;
00417 }
00418
00419
00420 theParticleChange.SetEnergyChange( sE );
00421 G4double mu = getMu( &anE_isoAng );
00422 theParticleChange.SetMomentumChange( 0.0 , std::sqrt ( 1 - mu*mu ) , mu );
00423
00424 }
00425
00426 else if ( random <= ( inelastic + theXSection->GetCoherentCrossSection( dp , theElement , theMaterial ) ) / total )
00427 {
00428
00429
00430 G4double E = aTrack.GetKineticEnergy();
00431
00432
00433 std::map < G4double , std::vector< std::pair< G4double , G4double >* >* >::iterator it;
00434 std::vector<G4double> v_temp;
00435 v_temp.clear();
00436 for ( it = coherentFSs.find( ielement )->second->begin() ; it != coherentFSs.find( ielement )->second->end() ; it++ )
00437 {
00438 v_temp.push_back( it->first );
00439 }
00440
00441
00442 std::pair < G4double , G4double > tempLH = find_LH ( aTemp , &v_temp );
00443
00444
00445
00446
00447 std::vector< std::pair< G4double , G4double >* >* pvE_p_TL = 0;
00448 std::vector< std::pair< G4double , G4double >* >* pvE_p_TH = 0;
00449
00450 if ( tempLH.first != 0.0 && tempLH.second != 0.0 )
00451 {
00452 pvE_p_TL = coherentFSs.find( ielement )->second->find ( tempLH.first/kelvin )->second;
00453 pvE_p_TH = coherentFSs.find( ielement )->second->find ( tempLH.first/kelvin )->second;
00454 }
00455 else if ( tempLH.first == 0.0 )
00456 {
00457 pvE_p_TL = coherentFSs.find( ielement )->second->find ( v_temp[ 0 ] )->second;
00458 pvE_p_TH = coherentFSs.find( ielement )->second->find ( v_temp[ 1 ] )->second;
00459 }
00460 else if ( tempLH.second == 0.0 )
00461 {
00462 pvE_p_TL = coherentFSs.find( ielement )->second->find ( v_temp.back() )->second;
00463 std::vector< G4double >::iterator itv;
00464 itv = v_temp.end();
00465 itv--;
00466 itv--;
00467 pvE_p_TL = coherentFSs.find( ielement )->second->find ( *itv )->second;
00468 }
00469
00470
00471 std::vector< G4double > vE_T;
00472 std::vector< G4double > vp_T;
00473
00474 G4int n1 = pvE_p_TL->size();
00475
00476
00477 for ( G4int i=1 ; i < n1 ; i++ )
00478 {
00479 if ( (*pvE_p_TL)[i]->first != (*pvE_p_TH)[i]->first ) G4HadronicException(__FILE__, __LINE__, "A problem is found in Thermal Scattering Data!");
00480 vE_T.push_back ( (*pvE_p_TL)[i]->first );
00481 vp_T.push_back ( get_linear_interpolated ( aTemp , std::pair< G4double , G4double > ( tempLH.first , (*pvE_p_TL)[i]->second ) , std::pair< G4double , G4double > ( tempLH.second , (*pvE_p_TL)[i]->second ) ) );
00482 }
00483
00484 G4int j = 0;
00485 for ( G4int i = 1 ; i < n ; i++ )
00486 {
00487 if ( E/eV < vE_T[ i ] )
00488 {
00489 j = i-1;
00490 break;
00491 }
00492 }
00493
00494 G4double rand_for_mu = G4UniformRand();
00495
00496 G4int k = 0;
00497 for ( G4int i = 1 ; i < j ; i++ )
00498 {
00499 G4double Pi = vp_T[ i ] / vp_T[ j ];
00500 if ( rand_for_mu < Pi )
00501 {
00502 k = i-1;
00503 break;
00504 }
00505 }
00506
00507
00508 G4double Ei = vE_T[ k ];
00509
00510 G4double mu = 1 - 2 * Ei / (E/eV) ;
00511
00512 if ( mu < -1.0 ) mu = -1.0;
00513
00514 theParticleChange.SetEnergyChange( E );
00515 theParticleChange.SetMomentumChange( 0.0 , std::sqrt ( 1 - mu*mu ) , mu );
00516
00517
00518 }
00519 else
00520 {
00521
00522
00523
00524 std::map < G4double , std::vector < E_isoAng* >* >::iterator it;
00525 std::vector<G4double> v_temp;
00526 v_temp.clear();
00527 for ( it = incoherentFSs.find( ielement )->second->begin() ; it != incoherentFSs.find( ielement )->second->end() ; it++ )
00528 {
00529 v_temp.push_back( it->first );
00530 }
00531
00532
00533 std::pair < G4double , G4double > tempLH = find_LH ( aTemp , &v_temp );
00534
00535
00536
00537
00538
00539 E_isoAng anEPM_TL_E;
00540 E_isoAng anEPM_TH_E;
00541
00542 if ( tempLH.first != 0.0 && tempLH.second != 0.0 )
00543 {
00544 anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( tempLH.first/kelvin )->second );
00545 anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( tempLH.second/kelvin )->second );
00546 }
00547 else if ( tempLH.first == 0.0 )
00548 {
00549 anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( v_temp[ 0 ] )->second );
00550 anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( v_temp[ 1 ] )->second );
00551 }
00552 else if ( tempLH.second == 0.0 )
00553 {
00554 anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( v_temp.back() )->second );
00555 std::vector< G4double >::iterator itv;
00556 itv = v_temp.end();
00557 itv--;
00558 itv--;
00559 anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( *itv )->second );
00560 }
00561
00562
00563 E_isoAng anEPM_T_E;
00564
00565 if ( anEPM_TL_E.n == anEPM_TH_E.n )
00566 {
00567 anEPM_T_E.n = anEPM_TL_E.n;
00568 for ( G4int i=0 ; i < anEPM_TL_E.n ; i++ )
00569 {
00570 G4double angle;
00571 angle = get_linear_interpolated ( aTemp , std::pair< G4double , G4double > ( tempLH.first , anEPM_TL_E.isoAngle[ i ] ) , std::pair< G4double , G4double > ( tempLH.second , anEPM_TH_E.isoAngle[ i ] ) );
00572 anEPM_T_E.isoAngle.push_back( angle );
00573 }
00574 }
00575 else
00576 {
00577 std::cout << "Do not Suuport yet." << std::endl;
00578 }
00579
00580
00581 G4double mu = getMu ( &anEPM_T_E );
00582
00583
00584 theParticleChange.SetEnergyChange( aTrack.GetKineticEnergy() );
00585 theParticleChange.SetMomentumChange( 0.0 , std::sqrt ( 1 - mu*mu ) , mu );
00586
00587 }
00588 delete dp;
00589
00590 return &theParticleChange;
00591
00592 }
00593 else
00594 {
00595
00596
00597 return theHPElastic -> ApplyYourself( aTrack, aNucleus );
00598 }
00599
00600 }
00601
00602
00603
00604 G4double G4NeutronHPThermalScattering::getMu( E_isoAng* anEPM )
00605 {
00606
00607 G4double random = G4UniformRand();
00608 G4double result = 0.0;
00609
00610 G4int in = int ( random * ( (*anEPM).n ) );
00611
00612 if ( in != 0 )
00613 {
00614 G4double mu_l = (*anEPM).isoAngle[ in-1 ];
00615 G4double mu_h = (*anEPM).isoAngle[ in ];
00616 result = ( mu_h - mu_l ) * ( random * ( (*anEPM).n ) - in ) + mu_l;
00617 }
00618 else
00619 {
00620 G4double x = random * (*anEPM).n;
00621 G4double D = ( (*anEPM).isoAngle[ 0 ] - ( -1 ) ) + ( 1 - (*anEPM).isoAngle[ (*anEPM).n - 1 ] );
00622 G4double ratio = ( (*anEPM).isoAngle[ 0 ] - ( -1 ) ) / D;
00623 if ( x <= ratio )
00624 {
00625 G4double mu_l = -1;
00626 G4double mu_h = (*anEPM).isoAngle[ 0 ];
00627 result = ( mu_h - mu_l ) * x + mu_l;
00628 }
00629 else
00630 {
00631 G4double mu_l = (*anEPM).isoAngle[ (*anEPM).n - 1 ];
00632 G4double mu_h = 1;
00633 result = ( mu_h - mu_l ) * x + mu_l;
00634 }
00635 }
00636 return result;
00637 }
00638
00639
00640
00641 std::pair < G4double , G4double > G4NeutronHPThermalScattering::find_LH ( G4double x , std::vector< G4double >* aVector )
00642 {
00643 G4double L = 0.0;
00644 G4double H = 0.0;
00645 std::vector< G4double >::iterator it;
00646 for ( it = aVector->begin() ; it != aVector->end() ; it++ )
00647 {
00648 if ( x <= *it )
00649 {
00650 H = *it;
00651 if ( it != aVector->begin() )
00652 {
00653 it--;
00654 L = *it;
00655 }
00656 else
00657 {
00658 L = 0.0;
00659 }
00660 break;
00661 }
00662 }
00663 if ( H == 0.0 )
00664 L = aVector->back();
00665
00666 return std::pair < G4double , G4double > ( L , H );
00667 }
00668
00669
00670
00671 G4double G4NeutronHPThermalScattering::get_linear_interpolated ( G4double x , std::pair< G4double , G4double > Low , std::pair< G4double , G4double > High )
00672 {
00673 G4double y=0.0;
00674 if ( High.first - Low.first != 0 )
00675 y = ( High.second - Low.second ) / ( High.first - Low.first ) * ( x - Low.first ) + Low.second;
00676 else
00677 std::cout << "G4NeutronHPThermalScattering liner interpolation err!!" << std::endl;
00678
00679 return y;
00680 }
00681
00682
00683
00684 E_isoAng G4NeutronHPThermalScattering::create_E_isoAng_from_energy ( G4double energy , std::vector< E_isoAng* >* vEPM )
00685 {
00686 E_isoAng anEPM_T_E;
00687
00688 std::vector< E_isoAng* >::iterator iv;
00689
00690 std::vector< G4double > v_e;
00691 v_e.clear();
00692 for ( iv = vEPM->begin() ; iv != vEPM->end() ; iv++ )
00693 v_e.push_back ( (*iv)->energy );
00694
00695 std::pair < G4double , G4double > energyLH = find_LH ( energy , &v_e );
00696
00697
00698 E_isoAng* panEPM_T_EL=0;
00699 E_isoAng* panEPM_T_EH=0;
00700
00701 if ( energyLH.first != 0.0 && energyLH.second != 0.0 )
00702 {
00703 for ( iv = vEPM->begin() ; iv != vEPM->end() ; iv++ )
00704 {
00705 if ( energyLH.first == (*iv)->energy )
00706 break;
00707 }
00708 panEPM_T_EL = *iv;
00709 iv++;
00710 panEPM_T_EH = *iv;
00711 }
00712 else if ( energyLH.first == 0.0 )
00713 {
00714 panEPM_T_EL = (*vEPM)[0];
00715 panEPM_T_EH = (*vEPM)[1];
00716 }
00717 else if ( energyLH.second == 0.0 )
00718 {
00719 panEPM_T_EH = (*vEPM).back();
00720 iv = vEPM->end();
00721 iv--;
00722 iv--;
00723 panEPM_T_EL = *iv;
00724 }
00725
00726 if ( panEPM_T_EL->n == panEPM_T_EH->n )
00727 {
00728 anEPM_T_E.energy = energy;
00729 anEPM_T_E.n = panEPM_T_EL->n;
00730
00731 for ( G4int i=0 ; i < panEPM_T_EL->n ; i++ )
00732 {
00733 G4double angle;
00734 angle = get_linear_interpolated ( energy , std::pair< G4double , G4double > ( energyLH.first , panEPM_T_EL->isoAngle[ i ] ) , std::pair< G4double , G4double > ( energyLH.second , panEPM_T_EH->isoAngle[ i ] ) );
00735 anEPM_T_E.isoAngle.push_back( angle );
00736 }
00737 }
00738 else
00739 {
00740 G4cout << "G4NeutronHPThermalScattering Do not Suuport yet." << G4endl;
00741 }
00742
00743
00744 return anEPM_T_E;
00745 }
00746
00747
00748
00749 G4double G4NeutronHPThermalScattering::get_secondary_energy_from_E_P_E_isoAng ( G4double random , E_P_E_isoAng* anE_P_E_isoAng )
00750 {
00751
00752 G4double secondary_energy = 0.0;
00753
00754 G4int n = anE_P_E_isoAng->n;
00755 G4double sum_p = 0.0;
00756 G4double sum_p_L = 0.0;
00757
00758 G4double total=0.0;
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772 total = anE_P_E_isoAng->sum_of_probXdEs;
00773
00774 for ( G4int i = 0 ; i < n-1 ; i++ )
00775 {
00776 G4double E_L = anE_P_E_isoAng->vE_isoAngle[i]->energy/eV;
00777 G4double E_H = anE_P_E_isoAng->vE_isoAngle[i+1]->energy/eV;
00778 G4double dE = E_H - E_L;
00779 sum_p += ( ( anE_P_E_isoAng->prob[i] ) * dE );
00780
00781 if ( random <= sum_p/total )
00782 {
00783 secondary_energy = get_linear_interpolated ( random , std::pair < G4double , G4double > ( sum_p_L/total , E_L ) , std::pair < G4double , G4double > ( sum_p/total , E_H ) );
00784 secondary_energy = secondary_energy*eV;
00785 break;
00786 }
00787 sum_p_L = sum_p;
00788 }
00789
00790 return secondary_energy;
00791 }
00792
00793
00794
00795 std::pair< G4double , E_isoAng > G4NeutronHPThermalScattering::create_sE_and_EPM_from_pE_and_vE_P_E_isoAng ( G4double rand_for_sE , G4double pE , std::vector < E_P_E_isoAng* >* vNEP_EPM )
00796 {
00797
00798 std::map< G4double , G4int > map_energy;
00799 map_energy.clear();
00800 std::vector< G4double > v_energy;
00801 v_energy.clear();
00802 std::vector< E_P_E_isoAng* >::iterator itv;
00803 G4int i = 0;
00804 for ( itv = vNEP_EPM->begin(); itv != vNEP_EPM->end(); itv++ )
00805 {
00806 v_energy.push_back( (*itv)->energy );
00807 map_energy.insert( std::pair < G4double , G4int > ( (*itv)->energy , i ) );
00808 i++;
00809 }
00810
00811 std::pair < G4double , G4double > energyLH = find_LH ( pE , &v_energy );
00812
00813 E_P_E_isoAng* pE_P_E_isoAng_EL = 0;
00814 E_P_E_isoAng* pE_P_E_isoAng_EH = 0;
00815
00816 if ( energyLH.first != 0.0 && energyLH.second != 0.0 )
00817 {
00818 pE_P_E_isoAng_EL = (*vNEP_EPM)[ map_energy.find ( energyLH.first )->second ];
00819 pE_P_E_isoAng_EH = (*vNEP_EPM)[ map_energy.find ( energyLH.second )->second ];
00820 }
00821 else if ( energyLH.first == 0.0 )
00822 {
00823 pE_P_E_isoAng_EL = (*vNEP_EPM)[ 0 ];
00824 pE_P_E_isoAng_EH = (*vNEP_EPM)[ 1 ];
00825 }
00826 if ( energyLH.second == 0.0 )
00827 {
00828 pE_P_E_isoAng_EH = (*vNEP_EPM).back();
00829 itv = vNEP_EPM->end();
00830 itv--;
00831 itv--;
00832 pE_P_E_isoAng_EL = *itv;
00833 }
00834
00835
00836 G4double sE;
00837 G4double sE_L;
00838 G4double sE_H;
00839
00840
00841 sE_L = get_secondary_energy_from_E_P_E_isoAng ( rand_for_sE , pE_P_E_isoAng_EL );
00842 sE_H = get_secondary_energy_from_E_P_E_isoAng ( rand_for_sE , pE_P_E_isoAng_EH );
00843
00844 sE = get_linear_interpolated ( pE , std::pair < G4double , G4double > ( energyLH.first , sE_L ) , std::pair < G4double , G4double > ( energyLH.second , sE_H ) );
00845
00846
00847 E_isoAng E_isoAng_L = create_E_isoAng_from_energy ( sE , &(pE_P_E_isoAng_EL->vE_isoAngle) );
00848 E_isoAng E_isoAng_H = create_E_isoAng_from_energy ( sE , &(pE_P_E_isoAng_EH->vE_isoAngle) );
00849
00850 E_isoAng anE_isoAng;
00851 if ( E_isoAng_L.n == E_isoAng_H.n )
00852 {
00853 anE_isoAng.n = E_isoAng_L.n;
00854 for ( G4int j=0 ; j < anE_isoAng.n ; j++ )
00855 {
00856 G4double angle;
00857 angle = get_linear_interpolated ( sE , std::pair< G4double , G4double > ( sE_L , E_isoAng_L.isoAngle[ j ] ) , std::pair< G4double , G4double > ( sE_H , E_isoAng_H.isoAngle[ j ] ) );
00858 anE_isoAng.isoAngle.push_back( angle );
00859 }
00860 }
00861 else
00862 {
00863 std::cout << "Do not Suuport yet." << std::endl;
00864 }
00865
00866
00867
00868 return std::pair< G4double , E_isoAng >( sE , anE_isoAng);
00869 }
00870
00871 void G4NeutronHPThermalScattering::buildPhysicsTable()
00872 {
00873
00874 dic.clear();
00875 std::map < G4String , G4int > co_dic;
00876
00877
00878 static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
00879 size_t numberOfMaterials = G4Material::GetNumberOfMaterials();
00880 for ( size_t i = 0 ; i < numberOfMaterials ; i++ )
00881 {
00882 G4Material* material = (*theMaterialTable)[i];
00883 size_t numberOfElements = material->GetNumberOfElements();
00884 for ( size_t j = 0 ; j < numberOfElements ; j++ )
00885 {
00886 const G4Element* element = material->GetElement(j);
00887 if ( names.IsThisThermalElement ( material->GetName() , element->GetName() ) )
00888 {
00889 G4int ts_ID_of_this_geometry;
00890 G4String ts_ndl_name = names.GetTS_NDL_Name( material->GetName() , element->GetName() );
00891 if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
00892 {
00893 ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) -> second;
00894 }
00895 else
00896 {
00897 ts_ID_of_this_geometry = co_dic.size();
00898 co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
00899 }
00900
00901
00902
00903
00904
00905 dic.insert( std::pair < std::pair < G4Material* , const G4Element* > , G4int > ( std::pair < G4Material* , const G4Element* > ( material , element ) , ts_ID_of_this_geometry ) );
00906 }
00907 }
00908 }
00909
00910
00911 static const G4ElementTable* theElementTable = G4Element::GetElementTable();
00912 size_t numberOfElements = G4Element::GetNumberOfElements();
00913
00914 for ( size_t i = 0 ; i < numberOfElements ; i++ )
00915 {
00916 const G4Element* element = (*theElementTable)[i];
00917 if ( names.IsThisThermalElement ( element->GetName() ) )
00918 {
00919 if ( names.IsThisThermalElement ( element->GetName() ) )
00920 {
00921 G4int ts_ID_of_this_geometry;
00922 G4String ts_ndl_name = names.GetTS_NDL_Name( element->GetName() );
00923 if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
00924 {
00925 ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) -> second;
00926 }
00927 else
00928 {
00929 ts_ID_of_this_geometry = co_dic.size();
00930 co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
00931 }
00932
00933
00934
00935
00936
00937 dic.insert( std::pair < std::pair < const G4Material* , const G4Element* > , G4int > ( std::pair < const G4Material* , const G4Element* > ( (G4Material*)NULL , element ) , ts_ID_of_this_geometry ) );
00938 }
00939 }
00940 }
00941
00942 G4cout << G4endl;
00943 G4cout << "Neutron HP Thermal Scattering: Following material-element pairs or elements are registered." << G4endl;
00944 for ( std::map < std::pair < const G4Material* , const G4Element* > , G4int >::iterator it = dic.begin() ; it != dic.end() ; it++ )
00945 {
00946 if ( it->first.first != NULL )
00947 {
00948 G4cout << "Material " << it->first.first->GetName() << " - Element " << it->first.second->GetName() << ", internal thermal scattering id " << it->second << G4endl;
00949 }
00950 else
00951 {
00952 G4cout << "Element " << it->first.second->GetName() << ", internal thermal scattering id " << it->second << G4endl;
00953 }
00954 }
00955 G4cout << G4endl;
00956
00957
00958
00959 G4String dirName;
00960 if ( !getenv( "G4NEUTRONHPDATA" ) )
00961 throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
00962 dirName = getenv( "G4NEUTRONHPDATA" );
00963
00964
00965
00966 G4String name;
00967
00968 for ( std::map < G4String , G4int >::iterator it = co_dic.begin() ; it != co_dic.end() ; it++ )
00969 {
00970 G4String tsndlName = it->first;
00971 G4int ts_ID = it->second;
00972
00973
00974 G4String fsName = "/ThermalScattering/Coherent/FS/";
00975 G4String fileName = dirName + fsName + tsndlName;
00976 coherentFSs.insert ( std::pair < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* > ( ts_ID , readACoherentFSDATA( fileName ) ) );
00977
00978
00979 fsName = "/ThermalScattering/Incoherent/FS/";
00980 fileName = dirName + fsName + tsndlName;
00981 incoherentFSs.insert ( std::pair < G4int , std::map < G4double , std::vector < E_isoAng* >* >* > ( ts_ID , readAnIncoherentFSDATA( fileName ) ) );
00982
00983
00984 fsName = "/ThermalScattering/Inelastic/FS/";
00985 fileName = dirName + fsName + tsndlName;
00986 inelasticFSs.insert ( std::pair < G4int , std::map < G4double , std::vector < E_P_E_isoAng* >* >* > ( ts_ID , readAnInelasticFSDATA( fileName ) ) );
00987 }
00988 }
00989
00990
00991 G4int G4NeutronHPThermalScattering::getTS_ID ( const G4Material* material , const G4Element* element )
00992 {
00993 G4int result = -1;
00994 if ( dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) ) != dic.end() )
00995 result = dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) )->second;
00996 return result;
00997 }
00998
00999 const std::pair<G4double, G4double> G4NeutronHPThermalScattering::GetFatalEnergyCheckLevels() const
01000 {
01001
01002 return std::pair<G4double, G4double>(10*perCent,DBL_MAX);
01003 }