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 #ifndef G4NeutronHPNBodyPhaseSpace_h
00030 #define G4NeutronHPNBodyPhaseSpace_h 1
00031
00032 #include <fstream>
00033 #include <CLHEP/Units/PhysicalConstants.h>
00034
00035 #include "globals.hh"
00036 #include "G4ios.hh"
00037 #include "G4Neutron.hh"
00038 #include "G4VNeutronHPEnergyAngular.hh"
00039 #include "G4ReactionProduct.hh"
00040
00041 class G4NeutronHPNBodyPhaseSpace : public G4VNeutronHPEnergyAngular
00042 {
00043 public:
00044
00045 G4NeutronHPNBodyPhaseSpace(){}
00046 ~G4NeutronHPNBodyPhaseSpace(){}
00047
00048 public:
00049
00050 void Init(G4double aMass, G4int aCount)
00051 {
00052 theTotalMass=aMass;
00053 theTotalCount=aCount;
00054 }
00055
00056 void Init(std::ifstream & aDataFile)
00057 {
00058 aDataFile >> theTotalMass >> theTotalCount;
00059 theTotalMass *= G4Neutron::Neutron()->GetPDGMass();
00060 }
00061
00062 G4ReactionProduct * Sample(G4double anEnergy, G4double massCode, G4double mass);
00063
00064 private:
00065
00066 inline G4double Prob(G4double anEnergy, G4double eMax, G4int n)
00067 {
00068 G4double result;
00069 result = std::sqrt(anEnergy)*std::pow(eMax-anEnergy, 3.*n/2.-4.);
00070 return result;
00071 }
00072
00073 inline G4double C(G4double anEnergy, G4double mass)
00074 {
00075 G4double result(0);
00076 if(theTotalCount==3) result = 4./CLHEP::pi/std::pow(GetEmax(anEnergy, mass),2);
00077 if(theTotalCount==4) result = 105./32./std::pow(GetEmax(anEnergy, mass), 3.5);
00078 if(theTotalCount==5) result = 256./14./CLHEP::pi/std::pow(GetEmax(anEnergy, mass), 5.);
00079 return result;
00080 }
00081
00082 inline G4double GetEmax(G4double anEnergy, G4double mass)
00083 {
00084 G4double result;
00085 G4double tMass = GetTarget()->GetMass();
00086 G4double pMass = GetNeutron()->GetMass();
00087 G4double availableEnergy = GetQValue() + anEnergy*tMass/(pMass+tMass);
00088 result = availableEnergy*(theTotalMass-mass)/theTotalMass;
00089 return result;
00090 }
00091
00092 G4double MeanEnergyOfThisInteraction() {return -1; }
00093
00094 private:
00095
00096 G4double theTotalMass;
00097 G4int theTotalCount;
00098
00099 };
00100 #endif