#include <G4LENDElastic.hh>
Inheritance diagram for G4LENDElastic:
Public Member Functions | |
G4LENDElastic (G4ParticleDefinition *pd) | |
~G4LENDElastic () | |
G4HadFinalState * | ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus) |
Definition at line 44 of file G4LENDElastic.hh.
G4LENDElastic::G4LENDElastic | ( | G4ParticleDefinition * | pd | ) | [inline] |
Definition at line 49 of file G4LENDElastic.hh.
References G4LENDModel::create_used_target_map(), and G4LENDModel::proj.
00050 :G4LENDModel( "LENDElastic" ) 00051 { 00052 proj = pd; 00053 00054 //theModelName = "LEND Elastic Model for "; 00055 //theModelName += proj->GetParticleName(); 00056 create_used_target_map(); 00057 };
G4LENDElastic::~G4LENDElastic | ( | ) | [inline] |
G4HadFinalState * G4LENDElastic::ApplyYourself | ( | const G4HadProjectile & | aTrack, | |
G4Nucleus & | aTargetNucleus | |||
) | [virtual] |
Reimplemented from G4LENDModel.
Definition at line 33 of file G4LENDElastic.cc.
References G4HadFinalState::AddSecondary(), G4HadFinalState::Clear(), G4ParticleTable::FindIon(), G4UniformRand, G4HadProjectile::Get4Momentum(), G4Nucleus::GetA_asInt(), G4HadProjectile::GetDefinition(), G4GIDI_target::getElasticFinalState(), G4ReactionProduct::GetKineticEnergy(), G4HadProjectile::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4HadProjectile::GetMaterial(), G4ReactionProduct::GetMomentum(), G4LENDManager::GetNucleusEncoding(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGMass(), G4Material::GetTemperature(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), G4LENDModel::lend_manager, G4ReactionProduct::Lorentz(), G4DynamicParticle::SetDefinition(), G4HadFinalState::SetEnergyChange(), G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMass(), G4DynamicParticle::SetMomentum(), G4ReactionProduct::SetMomentum(), G4HadFinalState::SetMomentumChange(), G4ReactionProduct::SetTotalEnergy(), G4HadronicInteraction::theParticleChange, and G4LENDModel::usedTarget_map.
00034 { 00035 00036 G4double temp = aTrack.GetMaterial()->GetTemperature(); 00037 00038 //G4int iZ = int ( aTarg.GetZ() ); 00039 //G4int iA = int ( aTarg.GetN() ); 00040 //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this) 00041 G4int iZ = aTarg.GetZ_asInt(); 00042 G4int iA = aTarg.GetA_asInt(); 00043 00044 G4double ke = aTrack.GetKineticEnergy(); 00045 00046 //G4HadFinalState* theResult = new G4HadFinalState(); 00047 G4HadFinalState* theResult = &theParticleChange; 00048 theResult->Clear(); 00049 00050 G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA ) )->second->GetTarget(); 00051 G4double aMu = aTarget->getElasticFinalState( ke*MeV, temp, NULL, NULL ); 00052 00053 G4double phi = twopi*G4UniformRand(); 00054 G4double theta = std::acos( aMu ); 00055 //G4double sinth = std::sin( theta ); 00056 00057 G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>( aTrack.GetDefinition() ) ); 00058 theNeutron.SetMomentum( aTrack.Get4Momentum().vect() ); 00059 theNeutron.SetKineticEnergy( ke ); 00060 00061 //G4cout << "iZ " << iZ << " iA " << iA << G4endl; 00062 00063 G4ReactionProduct theTarget( G4ParticleTable::GetParticleTable()->FindIon( iZ , iA , 0 , iZ ) ); 00064 00065 G4double mass = G4ParticleTable::GetParticleTable()->FindIon( iZ , iA , 0 , iZ )->GetPDGMass(); 00066 00067 // add Thermal motion 00068 G4double kT = k_Boltzmann*temp; 00069 G4ThreeVector v ( G4RandGauss::shoot() * std::sqrt( kT*mass ) 00070 , G4RandGauss::shoot() * std::sqrt( kT*mass ) 00071 , G4RandGauss::shoot() * std::sqrt( kT*mass ) ); 00072 theTarget.SetMomentum( v ); 00073 00074 G4ThreeVector the3Neutron = theNeutron.GetMomentum(); 00075 G4double nEnergy = theNeutron.GetTotalEnergy(); 00076 G4ThreeVector the3Target = theTarget.GetMomentum(); 00077 G4double tEnergy = theTarget.GetTotalEnergy(); 00078 G4ReactionProduct theCMS; 00079 G4double totE = nEnergy+tEnergy; 00080 G4ThreeVector the3CMS = the3Target+the3Neutron; 00081 theCMS.SetMomentum(the3CMS); 00082 G4double cmsMom = std::sqrt(the3CMS*the3CMS); 00083 G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom)); 00084 theCMS.SetMass(sqrts); 00085 theCMS.SetTotalEnergy(totE); 00086 00087 theNeutron.Lorentz(theNeutron, theCMS); 00088 theTarget.Lorentz(theTarget, theCMS); 00089 G4double en = theNeutron.GetTotalMomentum(); // already in CMS. 00090 G4ThreeVector cms3Mom=theNeutron.GetMomentum(); // for neutron direction in CMS 00091 G4double cms_theta=cms3Mom.theta(); 00092 G4double cms_phi=cms3Mom.phi(); 00093 G4ThreeVector tempVector; 00094 tempVector.setX( std::cos(theta)*std::sin(cms_theta)*std::cos(cms_phi) 00095 +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::cos(cms_phi) 00096 -std::sin(theta)*std::sin(phi)*std::sin(cms_phi) ); 00097 tempVector.setY( std::cos(theta)*std::sin(cms_theta)*std::sin(cms_phi) 00098 +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::sin(cms_phi) 00099 +std::sin(theta)*std::sin(phi)*std::cos(cms_phi) ); 00100 tempVector.setZ( std::cos(theta)*std::cos(cms_theta) 00101 -std::sin(theta)*std::cos(phi)*std::sin(cms_theta) ); 00102 tempVector *= en; 00103 theNeutron.SetMomentum(tempVector); 00104 theTarget.SetMomentum(-tempVector); 00105 G4double tP = theTarget.GetTotalMomentum(); 00106 G4double tM = theTarget.GetMass(); 00107 theTarget.SetTotalEnergy(std::sqrt((tP+tM)*(tP+tM)-2.*tP*tM)); 00108 00109 00110 theNeutron.Lorentz(theNeutron, -1.*theCMS); 00111 00112 //110913 Add Protection for very low energy (1e-6eV) scattering 00113 if ( theNeutron.GetKineticEnergy() <= 0 ) 00114 { 00115 theNeutron.SetTotalEnergy ( theNeutron.GetMass() * ( 1 + std::pow( 10 , -15.65 ) ) ); 00116 } 00117 00118 theTarget.Lorentz(theTarget, -1.*theCMS); 00119 if ( theTarget.GetKineticEnergy() < 0 ) 00120 { 00121 theTarget.SetTotalEnergy ( theTarget.GetMass() * ( 1 + std::pow( 10 , -15.65 ) ) ); 00122 } 00123 //110913 END 00124 00125 theTarget.Lorentz(theTarget, -1.*theCMS); 00126 00127 theResult->SetEnergyChange(theNeutron.GetKineticEnergy()); 00128 theResult->SetMomentumChange(theNeutron.GetMomentum().unit()); 00129 G4DynamicParticle* theRecoil = new G4DynamicParticle; 00130 00131 // theRecoil->SetDefinition( ionTable->GetIon( iZ , iA ) ); 00132 theRecoil->SetDefinition( G4ParticleTable::GetParticleTable()->FindIon( iZ, iA , 0, iZ )); 00133 theRecoil->SetMomentum( theTarget.GetMomentum() ); 00134 00135 theResult->AddSecondary( theRecoil ); 00136 00137 return theResult; 00138 00139 }