Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4ParticleHPThermalScattering Class Reference

#include <G4ParticleHPThermalScattering.hh>

+ Inheritance diagram for G4ParticleHPThermalScattering:

Public Member Functions

 G4ParticleHPThermalScattering ()
 
 ~G4ParticleHPThermalScattering ()
 
G4HadFinalState * ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus)
 
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels () const
 
void AddUserThermalScatteringFile (G4String, G4String)
 
void BuildPhysicsTable (const G4ParticleDefinition &)
 
virtual void ModelDescription (std::ostream &outFile) const
 
- Public Member Functions inherited from G4HadronicInteraction
 G4HadronicInteraction (const G4String &modelName="HadronicModel")
 
virtual ~G4HadronicInteraction ()
 
virtual G4HadFinalState * ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
virtual G4bool IsApplicable (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
G4double GetMinEnergy () const
 
G4double GetMinEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMinEnergy (G4double anEnergy)
 
void SetMinEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMinEnergy (G4double anEnergy, const G4Material *aMaterial)
 
G4double GetMaxEnergy () const
 
G4double GetMaxEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMaxEnergy (const G4double anEnergy)
 
void SetMaxEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMaxEnergy (G4double anEnergy, const G4Material *aMaterial)
 
G4int GetVerboseLevel () const
 
void SetVerboseLevel (G4int value)
 
const G4String & GetModelName () const
 
void DeActivateFor (const G4Material *aMaterial)
 
void ActivateFor (const G4Material *aMaterial)
 
void DeActivateFor (const G4Element *anElement)
 
void ActivateFor (const G4Element *anElement)
 
G4bool IsBlocked (const G4Material *aMaterial) const
 
G4bool IsBlocked (const G4Element *anElement) const
 
void SetRecoilEnergyThreshold (G4double val)
 
G4double GetRecoilEnergyThreshold () const
 
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels () const
 
virtual std::pair< G4double, G4double > GetEnergyMomentumCheckLevels () const
 
void SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
 
virtual void ModelDescription (std::ostream &outFile) const
 
virtual void BuildPhysicsTable (const G4ParticleDefinition &)
 
virtual void InitialiseModel ()
 
 G4HadronicInteraction (const G4HadronicInteraction &right)=delete
 
const G4HadronicInteraction & operator= (const G4HadronicInteraction &right)=delete
 
G4bool operator== (const G4HadronicInteraction &right) const =delete
 
G4bool operator!= (const G4HadronicInteraction &right) const =delete
 

Additional Inherited Members

- Protected Member Functions inherited from G4HadronicInteraction
void SetModelName (const G4String &nam)
 
G4bool IsBlocked () const
 
void Block ()
 
- Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState theParticleChange
 
G4int verboseLevel
 
G4double theMinEnergy
 
G4double theMaxEnergy
 
G4bool isBlocked
 

Detailed Description

Definition at line 82 of file G4ParticleHPThermalScattering.hh.

Constructor & Destructor Documentation

◆ G4ParticleHPThermalScattering()

G4ParticleHPThermalScattering::G4ParticleHPThermalScattering ( )

Definition at line 56 of file G4ParticleHPThermalScattering.cc.

57 :G4HadronicInteraction("NeutronHPThermalScattering")
58,coherentFSs(nullptr)
59,incoherentFSs(nullptr)
60,inelasticFSs(nullptr)
61{
62 theHPElastic = new G4ParticleHPElastic();
63
64 SetMinEnergy( 0.*eV );
65 SetMaxEnergy( 4*eV );
66 theXSection = new G4ParticleHPThermalScatteringData();
67
68 nMaterial = 0;
69 nElement = 0;
70}
G4HadronicInteraction::SetMinEnergy
void SetMinEnergy(G4double anEnergy)
Definition: G4HadronicInteraction.hh:89
G4HadronicInteraction::SetMaxEnergy
void SetMaxEnergy(const G4double anEnergy)
Definition: G4HadronicInteraction.hh:102

◆ ~G4ParticleHPThermalScattering()

G4ParticleHPThermalScattering::~G4ParticleHPThermalScattering ( )

Definition at line 73 of file G4ParticleHPThermalScattering.cc.

74{
75 delete theHPElastic;
76}

Member Function Documentation

◆ AddUserThermalScatteringFile()

void G4ParticleHPThermalScattering::AddUserThermalScatteringFile ( G4String  nameG4Element,
G4String  filename 
)

Definition at line 1227 of file G4ParticleHPThermalScattering.cc.

1228{
1229 names.AddThermalElement( nameG4Element , filename );
1230 theXSection->AddUserThermalScatteringFile( nameG4Element , filename );
1231 buildPhysicsTable();
1232}

◆ ApplyYourself()

G4HadFinalState * G4ParticleHPThermalScattering::ApplyYourself ( const G4HadProjectile &  aTrack,
G4Nucleus &  aTargetNucleus 
)
virtual

Reimplemented from G4HadronicInteraction.

Definition at line 335 of file G4ParticleHPThermalScattering.cc.

336{
337
338// Select Element > Reaction >
339
340 const G4Material * theMaterial = aTrack.GetMaterial();
341 G4double aTemp = theMaterial->GetTemperature();
342 G4int n = (G4int)theMaterial->GetNumberOfElements();
343
344 G4bool findThermalElement = false;
345 G4int ielement;
346 const G4Element* theElement = nullptr;
347 for ( G4int i = 0; i < n ; ++i )
348 {
349 theElement = theMaterial->GetElement(i);
350 // Select target element
351 if ( aNucleus.GetZ_asInt() == (G4int)(theElement->GetZ() + 0.5 ) )
352 {
353 //Check Applicability of Thermal Scattering
354 if ( getTS_ID( nullptr , theElement ) != -1 )
355 {
356 ielement = getTS_ID( nullptr , theElement );
357 findThermalElement = true;
358 break;
359 }
360 else if ( getTS_ID( theMaterial , theElement ) != -1 )
361 {
362 ielement = getTS_ID( theMaterial , theElement );
363 findThermalElement = true;
364 break;
365 }
366 }
367 }
368
369 if ( findThermalElement == true )
370 {
371
372 // Select Reaction (Inelastic, coherent, incoherent)
373 const G4ParticleDefinition* pd = aTrack.GetDefinition();
374 G4DynamicParticle* dp = new G4DynamicParticle ( pd , aTrack.Get4Momentum() );
375 G4double total = theXSection->GetCrossSection( dp , theElement , theMaterial );
376 G4double inelastic = theXSection->GetInelasticCrossSection( dp , theElement , theMaterial );
377
378 G4double random = G4UniformRand();
379 if ( random <= inelastic/total )
380 {
381 // Inelastic
382
383 std::vector<G4double> v_temp;
384 v_temp.clear();
385 for (auto it = inelasticFSs->find( ielement )->second->cbegin();
386 it != inelasticFSs->find( ielement )->second->cend() ; ++it )
387 {
388 v_temp.push_back( it->first );
389 }
390
391 std::pair < G4double , G4double > tempLH = find_LH ( aTemp , &v_temp );
392 //
393 // For T_L aNEP_EPM_TL and T_H aNEP_EPM_TH
394 //
395 std::vector< E_P_E_isoAng* >* vNEP_EPM_TL = nullptr;
396 std::vector< E_P_E_isoAng* >* vNEP_EPM_TH = nullptr;
397
398 if ( tempLH.first != 0.0 && tempLH.second != 0.0 )
399 {
400 vNEP_EPM_TL = inelasticFSs->find( ielement )->second->find ( tempLH.first/kelvin )->second;
401 vNEP_EPM_TH = inelasticFSs->find( ielement )->second->find ( tempLH.second/kelvin )->second;
402 }
403 else if ( tempLH.first == 0.0 )
404 {
405 auto itm = inelasticFSs->find( ielement )->second->cbegin();
406 vNEP_EPM_TL = itm->second;
407 ++itm;
408 vNEP_EPM_TH = itm->second;
409 tempLH.first = tempLH.second;
410 tempLH.second = itm->first;
411 }
412 else if ( tempLH.second == 0.0 )
413 {
414 auto itm = inelasticFSs->find( ielement )->second->cend();
415 --itm;
416 vNEP_EPM_TH = itm->second;
417 --itm;
418 vNEP_EPM_TL = itm->second;
419 tempLH.second = tempLH.first;
420 tempLH.first = itm->first;
421 }
422
423 G4double sE=0., mu=1.0;
424
425 // New Geant4 method - Stochastic temperature interpolation of the final state
426 // (continuous temperature interpolation was used previously)
427 std::pair< G4double , G4double > secondaryParam;
428 G4double rand_temp = G4UniformRand();
429 if ( rand_temp < (aTemp-tempLH.first)/(tempLH.second - tempLH.first) )
430 secondaryParam = sample_inelastic_E_mu( aTrack.GetKineticEnergy() , vNEP_EPM_TH );
431 else
432 secondaryParam = sample_inelastic_E_mu( aTrack.GetKineticEnergy() , vNEP_EPM_TL );
433
434 sE = secondaryParam.first;
435 mu = secondaryParam.second;
436
437 //set
438 theParticleChange.SetEnergyChange( sE );
439 G4double phi = CLHEP::twopi*G4UniformRand();
440 G4double sint= std::sqrt ( 1 - mu*mu );
441 theParticleChange.SetMomentumChange( sint*std::cos(phi), sint*std::sin(phi), mu );
442 }
443 else if ( random <= ( inelastic + theXSection->GetCoherentCrossSection( dp , theElement , theMaterial ) ) / total )
444 {
445 // Coherent Elastic
446
447 G4double E = aTrack.GetKineticEnergy();
448
449 // T_L and T_H
450 std::vector<G4double> v_temp;
451 v_temp.clear();
452 for (auto it = coherentFSs->find(ielement)->second->cbegin();
453 it != coherentFSs->find(ielement)->second->cend(); ++it)
454 {
455 v_temp.push_back( it->first );
456 }
457
458 // T_L T_H
459 std::pair < G4double , G4double > tempLH = find_LH ( aTemp , &v_temp );
460 //
461 //
462 // For T_L anEPM_TL and T_H anEPM_TH
463 //
464 std::vector< std::pair< G4double , G4double >* >* pvE_p_TL = nullptr;
465 std::vector< std::pair< G4double , G4double >* >* pvE_p_TH = nullptr;
466
467 if ( tempLH.first != 0.0 && tempLH.second != 0.0 )
468 {
469 pvE_p_TL = coherentFSs->find( ielement )->second->find ( tempLH.first/kelvin )->second;
470 pvE_p_TH = coherentFSs->find( ielement )->second->find ( tempLH.first/kelvin )->second;
471 }
472 else if ( tempLH.first == 0.0 )
473 {
474 pvE_p_TL = coherentFSs->find( ielement )->second->find ( v_temp[ 0 ] )->second;
475 pvE_p_TH = coherentFSs->find( ielement )->second->find ( v_temp[ 1 ] )->second;
476 tempLH.first = tempLH.second;
477 tempLH.second = v_temp[ 1 ];
478 }
479 else if ( tempLH.second == 0.0 )
480 {
481 pvE_p_TH = coherentFSs->find( ielement )->second->find ( v_temp.back() )->second;
482 auto itv = v_temp.cend();
483 --itv;
484 --itv;
485 pvE_p_TL = coherentFSs->find( ielement )->second->find ( *itv )->second;
486 tempLH.second = tempLH.first;
487 tempLH.first = *itv;
488 }
489 else
490 {
491 // tempLH.first == 0.0 && tempLH.second
492 throw G4HadronicException(__FILE__, __LINE__, "A problem is found in Thermal Scattering Data! Unexpected temperature values in data");
493 }
494
495 std::vector< G4double > vE_T;
496 std::vector< G4double > vp_T;
497
498 G4int n1 = (G4int)pvE_p_TL->size();
499
500 // New Geant4 method - Stochastic interpolation of the final state
501 std::vector< std::pair< G4double , G4double >* >* pvE_p_T_sampled;
502 G4double rand_temp = G4UniformRand();
503 if ( rand_temp < (aTemp-tempLH.first)/(tempLH.second - tempLH.first) )
504 pvE_p_T_sampled = pvE_p_TH;
505 else
506 pvE_p_T_sampled = pvE_p_TL;
507
508 //171005 fix bug, contribution from H.N. TRAN@CEA
509 for ( G4int i=0 ; i < n1 ; ++i )
510 {
511 vE_T.push_back ( (*pvE_p_T_sampled)[i]->first );
512 vp_T.push_back ( (*pvE_p_T_sampled)[i]->second );
513 }
514
515 G4int j = 0;
516 for ( G4int i = 1 ; i < n1 ; ++i )
517 {
518 if ( E/eV < vE_T[ i ] )
519 {
520 j = i-1;
521 break;
522 }
523 }
524
525 G4double rand_for_mu = G4UniformRand();
526
527 G4int k = 0;
528 for ( G4int i = 0 ; i <= j ; ++i )
529 {
530 G4double Pi = vp_T[ i ] / vp_T[ j ];
531 if ( rand_for_mu < Pi )
532 {
533 k = i;
534 break;
535 }
536 }
537
538 G4double Ei = vE_T[ k ];
539
540 G4double mu = 1 - 2 * Ei / (E/eV) ;
541
542 if ( mu < -1.0 ) mu = -1.0;
543
544 theParticleChange.SetEnergyChange( E );
545 G4double phi = CLHEP::twopi*G4UniformRand();
546 G4double sint= std::sqrt ( 1 - mu*mu );
547 theParticleChange.SetMomentumChange( sint*std::cos(phi), sint*std::sin(phi), mu );
548 }
549 else
550 {
551 // InCoherent Elastic
552
553 // T_L and T_H
554 std::vector<G4double> v_temp;
555 v_temp.clear();
556 for (auto it = incoherentFSs->find(ielement)->second->cbegin();
557 it != incoherentFSs->find(ielement)->second->cend(); ++it)
558 {
559 v_temp.push_back( it->first );
560 }
561
562 // T_L T_H
563 std::pair < G4double , G4double > tempLH = find_LH ( aTemp , &v_temp );
564
565 //
566 // For T_L anEPM_TL and T_H anEPM_TH
567 //
568
569 E_isoAng anEPM_TL_E;
570 E_isoAng anEPM_TH_E;
571
572 if ( tempLH.first != 0.0 && tempLH.second != 0.0 ) {
573 //Interpolate TL and TH
574 anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs->find( ielement )->second->find ( tempLH.first/kelvin )->second );
575 anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs->find( ielement )->second->find ( tempLH.second/kelvin )->second );
576 } else if ( tempLH.first == 0.0 ) {
577 //Extrapolate T0 and T1
578 anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs->find( ielement )->second->find ( v_temp[ 0 ] )->second );
579 anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs->find( ielement )->second->find ( v_temp[ 1 ] )->second );
580 tempLH.first = tempLH.second;
581 tempLH.second = v_temp[ 1 ];
582 } else if ( tempLH.second == 0.0 ) {
583 //Extrapolate Tmax-1 and Tmax
584 anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs->find( ielement )->second->find ( v_temp.back() )->second );
585 auto itv = v_temp.cend();
586 --itv;
587 --itv;
588 anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs->find( ielement )->second->find ( *itv )->second );
589 tempLH.second = tempLH.first;
590 tempLH.first = *itv;
591 }
592
593 // E_isoAng for aTemp and aTrack.GetKineticEnergy()
594 G4double mu=1.0;
595
596 // New Geant4 method - Stochastic interpolation of the final state
597 E_isoAng anEPM_T_E_sampled;
598 G4double rand_temp = G4UniformRand();
599 if ( rand_temp < (aTemp-tempLH.first)/(tempLH.second - tempLH.first) )
600 anEPM_T_E_sampled = anEPM_TH_E;
601 else
602 anEPM_T_E_sampled = anEPM_TL_E;
603
604 mu = getMu ( &anEPM_T_E_sampled );
605
606 // Set Final State
607 theParticleChange.SetEnergyChange( aTrack.GetKineticEnergy() ); // No energy change in Elastic
608 G4double phi = CLHEP::twopi*G4UniformRand();
609 G4double sint= std::sqrt ( 1 - mu*mu );
610 theParticleChange.SetMomentumChange( sint*std::cos(phi), sint*std::sin(phi), mu );
611 }
612 delete dp;
613
614 return &theParticleChange;
615 }
616 else
617 {
618 // Not thermal element
619 // Neutron HP will handle
620 return theHPElastic -> ApplyYourself( aTrack, aNucleus, 1); // L. Thulliez 2021/05/04 (CEA-Saclay)
621 }
622}
G4double
double G4double
Definition: G4Types.hh:83
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
G4Element::GetZ
G4double GetZ() const
Definition: G4Element.hh:131
G4HadFinalState::SetEnergyChange
void SetEnergyChange(G4double anEnergy)
Definition: G4HadFinalState.cc:39
G4HadFinalState::SetMomentumChange
void SetMomentumChange(const G4ThreeVector &aV)
Definition: G4HadFinalState.hh:56
G4HadProjectile::GetMaterial
const G4Material * GetMaterial() const
Definition: G4HadProjectile.hh:81
G4HadProjectile::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4HadProjectile.hh:86
G4HadProjectile::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4HadProjectile.hh:116
G4HadProjectile::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4HadProjectile.hh:91
G4Material::GetTemperature
G4double GetTemperature() const
Definition: G4Material.hh:177
G4Material::GetElement
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:197
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:181
G4ParticleHPThermalScatteringData::GetCrossSection
G4double GetCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
Definition: G4ParticleHPThermalScatteringData.cc:368
G4ParticleHPThermalScatteringData::GetCoherentCrossSection
G4double GetCoherentCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
Definition: G4ParticleHPThermalScatteringData.cc:397
G4ParticleHPThermalScatteringData::GetInelasticCrossSection
G4double GetInelasticCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
Definition: G4ParticleHPThermalScatteringData.cc:388
G4ParticleHPThermalScattering::ApplyYourself
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus)
Definition: G4ParticleHPThermalScattering.cc:335
G4INCL::CrossSections::total
G4double total(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:60
E_isoAng
Definition: G4ParticleHPThermalScattering.hh:53

Referenced by ApplyYourself().

◆ BuildPhysicsTable()

void G4ParticleHPThermalScattering::BuildPhysicsTable ( const G4ParticleDefinition &  particle)
virtual

Reimplemented from G4HadronicInteraction.

Definition at line 139 of file G4ParticleHPThermalScattering.cc.

140{
141 buildPhysicsTable();
142 theHPElastic->BuildPhysicsTable( particle );
143}
G4ParticleHPElastic::BuildPhysicsTable
void BuildPhysicsTable(const G4ParticleDefinition &)
Definition: G4ParticleHPElastic.cc:157

◆ GetFatalEnergyCheckLevels()

const std::pair< G4double, G4double > G4ParticleHPThermalScattering::GetFatalEnergyCheckLevels ( ) const
virtual

Reimplemented from G4HadronicInteraction.

Definition at line 1220 of file G4ParticleHPThermalScattering.cc.

1221{
1222 // max energy non-conservation is mass of heavy nucleus
1223 return std::pair<G4double, G4double>(10.0*perCent, 350.0*CLHEP::GeV);
1224}

◆ ModelDescription()

void G4ParticleHPThermalScattering::ModelDescription ( std::ostream &  outFile) const
virtual

Reimplemented from G4HadronicInteraction.

Definition at line 1250 of file G4ParticleHPThermalScattering.cc.

1251{
1252 outFile << "High Precision model based on thermal scattering data in\n"
1253 << "evaluated nuclear data libraries for neutrons below 5eV\n"
1254 << "on specific materials\n";
1255}
The documentation for this class was generated from the following files: