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 G4NeutronHPMadlandNixSpectrum_h
00030 #define G4NeutronHPMadlandNixSpectrum_h 1
00031
00032 #include <fstream>
00033 #include <cmath>
00034 #include <CLHEP/Units/PhysicalConstants.h>
00035
00036 #include "globals.hh"
00037 #include "G4ios.hh"
00038 #include "Randomize.hh"
00039 #include "G4NeutronHPVector.hh"
00040 #include "G4VNeutronHPEDis.hh"
00041
00042
00043
00044
00045
00046
00047
00048 class G4NeutronHPMadlandNixSpectrum : public G4VNeutronHPEDis
00049 {
00050 public:
00051 G4NeutronHPMadlandNixSpectrum()
00052 {
00053 expm1 = std::exp(-1.);
00054 }
00055 ~G4NeutronHPMadlandNixSpectrum()
00056 {
00057 }
00058
00059 inline void Init(std::ifstream & aDataFile)
00060 {
00061 theFractionalProb.Init(aDataFile);
00062 aDataFile>> theAvarageKineticPerNucleonForLightFragments;
00063 theAvarageKineticPerNucleonForLightFragments*=CLHEP::eV;
00064 aDataFile>> theAvarageKineticPerNucleonForHeavyFragments;
00065 theAvarageKineticPerNucleonForHeavyFragments*=CLHEP::eV;
00066 theMaxTemp.Init(aDataFile);
00067 }
00068
00069 inline G4double GetFractionalProbability(G4double anEnergy)
00070 {
00071 return theFractionalProb.GetY(anEnergy);
00072 }
00073
00074 G4double Sample(G4double anEnergy);
00075
00076 private:
00077
00078 G4double Madland(G4double aSecEnergy, G4double tm);
00079
00080 inline G4double FissionIntegral(G4double tm, G4double anEnergy)
00081 {
00082 return 0.5*( GIntegral(tm, anEnergy, theAvarageKineticPerNucleonForLightFragments)
00083 +GIntegral(tm, anEnergy, theAvarageKineticPerNucleonForHeavyFragments) );
00084 }
00085
00086 G4double GIntegral(G4double tm, G4double anEnergy, G4double aMean);
00087
00088 inline G4double Gamma05(G4double aValue)
00089 {
00090 G4double result;
00091
00092 G4double x = std::sqrt(aValue);
00093 G4double t = 1./(1+0.47047*x);
00094 result = 1- (0.3480242*t - 0.0958798*t*t + 0.7478556*t*t*t)*std::exp(-aValue);
00095 result *= std::sqrt(CLHEP::pi);
00096 return result;
00097 }
00098
00099 inline G4double Gamma15(G4double aValue)
00100 {
00101 G4double result;
00102
00103 result = 0.5*Gamma05(aValue) - std::sqrt(aValue)*std::exp(-aValue);
00104 return result;
00105 }
00106
00107 inline G4double Gamma25(G4double aValue)
00108 {
00109 G4double result;
00110 result = 1.5*Gamma15(aValue) - std::pow(aValue,1.5)*std::exp(aValue);
00111 return result;
00112 }
00113
00114 inline G4double E1(G4double aValue)
00115 {
00116
00117
00118 G4double gamma = 0.577216;
00119 G4double precision = 0.000001;
00120 G4double result =-gamma - std::log(aValue);
00121 G4double term = -aValue;
00122
00123
00124 G4int count = 1;
00125 result -= term;
00126 for(;;)
00127 {
00128 count++;
00129
00130
00131 term = -term*aValue*(count-1)/(count*count);
00132 result -=term;
00133 if(std::fabs(term)/std::fabs(result)<precision) break;
00134 }
00135
00136
00137 return result;
00138 }
00139
00140 private:
00141
00142 G4double expm1;
00143
00144 private:
00145
00146 G4NeutronHPVector theFractionalProb;
00147
00148 G4double theAvarageKineticPerNucleonForLightFragments;
00149 G4double theAvarageKineticPerNucleonForHeavyFragments;
00150
00151 G4NeutronHPVector theMaxTemp;
00152
00153 };
00154
00155 #endif