#include <G4NeutronHPNBodyPhaseSpace.hh>
Inheritance diagram for G4NeutronHPNBodyPhaseSpace:
Public Member Functions | |
G4NeutronHPNBodyPhaseSpace () | |
~G4NeutronHPNBodyPhaseSpace () | |
void | Init (G4double aMass, G4int aCount) |
void | Init (std::ifstream &aDataFile) |
G4ReactionProduct * | Sample (G4double anEnergy, G4double massCode, G4double mass) |
Definition at line 41 of file G4NeutronHPNBodyPhaseSpace.hh.
G4NeutronHPNBodyPhaseSpace::G4NeutronHPNBodyPhaseSpace | ( | ) | [inline] |
G4NeutronHPNBodyPhaseSpace::~G4NeutronHPNBodyPhaseSpace | ( | ) | [inline] |
void G4NeutronHPNBodyPhaseSpace::Init | ( | std::ifstream & | aDataFile | ) | [inline, virtual] |
Implements G4VNeutronHPEnergyAngular.
Definition at line 56 of file G4NeutronHPNBodyPhaseSpace.hh.
References G4ParticleDefinition::GetPDGMass(), and G4Neutron::Neutron().
00057 { 00058 aDataFile >> theTotalMass >> theTotalCount; 00059 theTotalMass *= G4Neutron::Neutron()->GetPDGMass(); 00060 }
Definition at line 50 of file G4NeutronHPNBodyPhaseSpace.hh.
Referenced by G4NeutronHPInelasticBaseFS::BaseApply().
G4ReactionProduct * G4NeutronHPNBodyPhaseSpace::Sample | ( | G4double | anEnergy, | |
G4double | massCode, | |||
G4double | mass | |||
) | [virtual] |
Implements G4VNeutronHPEnergyAngular.
Definition at line 43 of file G4NeutronHPNBodyPhaseSpace.cc.
References G4Alpha::Alpha(), G4Deuteron::Deuteron(), G4Electron::Electron(), G4UniformRand, G4Gamma::Gamma(), G4ReactionProduct::GetMass(), G4VNeutronHPEnergyAngular::GetNeutron(), G4VNeutronHPEnergyAngular::GetTarget(), G4ReactionProduct::GetTotalMomentum(), G4He3::He3(), G4ReactionProduct::Lorentz(), G4Neutron::Neutron(), G4Positron::Positron(), G4Proton::Proton(), G4ReactionProduct::SetDefinition(), G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMomentum(), and G4Triton::Triton().
Referenced by G4NeutronHPInelasticBaseFS::BaseApply().
00044 { 00045 G4ReactionProduct * result = new G4ReactionProduct; 00046 G4int Z = static_cast<G4int>(massCode/1000); 00047 G4int A = static_cast<G4int>(massCode-1000*Z); 00048 00049 if(massCode==0) 00050 { 00051 result->SetDefinition(G4Gamma::Gamma()); 00052 } 00053 else if(A==0) 00054 { 00055 result->SetDefinition(G4Electron::Electron()); 00056 if(Z==1) result->SetDefinition(G4Positron::Positron()); 00057 } 00058 else if(A==1) 00059 { 00060 result->SetDefinition(G4Neutron::Neutron()); 00061 if(Z==1) result->SetDefinition(G4Proton::Proton()); 00062 } 00063 else if(A==2) 00064 { 00065 result->SetDefinition(G4Deuteron::Deuteron()); 00066 } 00067 else if(A==3) 00068 { 00069 result->SetDefinition(G4Triton::Triton()); 00070 if(Z==2) result->SetDefinition(G4He3::He3()); 00071 } 00072 else if(A==4) 00073 { 00074 result->SetDefinition(G4Alpha::Alpha()); 00075 if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "Unknown ion case 1"); 00076 } 00077 else 00078 { 00079 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPNBodyPhaseSpace: Unknown ion case 2"); 00080 } 00081 00082 // Get the energy from phase-space distribution 00083 // in CMS 00084 // P = Cn*std::sqrt(E')*(Emax-E')**(3*n/2-4) 00085 G4double maxE = GetEmax(anEnergy, result->GetMass()); 00086 G4double energy; 00087 G4double max(0); 00088 if(theTotalCount<=3) 00089 { 00090 max = maxE/2.; 00091 } 00092 else if(theTotalCount==4) 00093 { 00094 max = maxE/5.; 00095 } 00096 else if(theTotalCount==5) 00097 { 00098 max = maxE/8.; 00099 } 00100 else 00101 { 00102 throw G4HadronicException(__FILE__, __LINE__, "NeutronHP Phase-space distribution cannot cope with this number of particles"); 00103 } 00104 G4double testit; 00105 G4double rand0 = Prob(max, maxE, theTotalCount); 00106 G4double rand; 00107 00108 do 00109 { 00110 rand = rand0*G4UniformRand(); 00111 energy = maxE*G4UniformRand(); 00112 testit = Prob(energy, maxE, theTotalCount); 00113 } 00114 while(rand > testit); 00115 result->SetKineticEnergy(energy); 00116 00117 // now do random direction 00118 G4double cosTh = 2.*G4UniformRand()-1.; 00119 G4double phi = twopi*G4UniformRand(); 00120 G4double theta = std::acos(cosTh); 00121 G4double sinth = std::sin(theta); 00122 G4double mtot = result->GetTotalMomentum(); 00123 G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) ); 00124 result->SetMomentum(tempVector); 00125 G4ReactionProduct aCMS = *GetTarget()+*GetNeutron(); 00126 result->Lorentz(*result, -1.*aCMS); 00127 return result; 00128 }