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 #include "G4NeutronHPNBodyPhaseSpace.hh"
00030 #include "G4PhysicalConstants.hh"
00031 #include "Randomize.hh"
00032 #include "G4ThreeVector.hh"
00033 #include "G4Gamma.hh"
00034 #include "G4Electron.hh"
00035 #include "G4Positron.hh"
00036 #include "G4Neutron.hh"
00037 #include "G4Proton.hh"
00038 #include "G4Deuteron.hh"
00039 #include "G4Triton.hh"
00040 #include "G4He3.hh"
00041 #include "G4Alpha.hh"
00042
00043 G4ReactionProduct * G4NeutronHPNBodyPhaseSpace::Sample(G4double anEnergy, G4double massCode, G4double )
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
00083
00084
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
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 }