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
00030 #include "G4NeutronHPMadlandNixSpectrum.hh"
00031 #include "G4SystemOfUnits.hh"
00032
00033 G4double G4NeutronHPMadlandNixSpectrum::Madland(G4double aSecEnergy, G4double tm)
00034 {
00035 G4double result;
00036 G4double energy = aSecEnergy/eV;
00037 G4double EF;
00038
00039 EF = theAvarageKineticPerNucleonForLightFragments/eV;
00040 G4double lightU1 = std::sqrt(energy)-std::sqrt(EF);
00041 lightU1 *= lightU1/tm;
00042 G4double lightU2 = std::sqrt(energy)+std::sqrt(EF);
00043 lightU2 *= lightU2/tm;
00044 G4double lightTerm=0;
00045 if(theAvarageKineticPerNucleonForLightFragments>1*eV)
00046 {
00047 lightTerm = std::pow(lightU2, 1.5)*E1(lightU2);
00048 lightTerm -= std::pow(lightU1, 1.5)*E1(lightU1);
00049 lightTerm += Gamma15(lightU2)-Gamma15(lightU1);
00050 lightTerm /= 3.*std::sqrt(tm*EF);
00051 }
00052
00053 EF = theAvarageKineticPerNucleonForHeavyFragments/eV;
00054 G4double heavyU1 = std::sqrt(energy)-std::sqrt(EF);
00055 heavyU1 *= heavyU1/tm;
00056 G4double heavyU2 = std::sqrt(energy)+std::sqrt(EF);
00057 heavyU2 *= heavyU2/tm;
00058 G4double heavyTerm=0 ;
00059 if(theAvarageKineticPerNucleonForHeavyFragments> 1*eV)
00060 {
00061 heavyTerm = std::pow(heavyU2, 1.5)*E1(heavyU2);
00062 heavyTerm -= std::pow(heavyU1, 1.5)*E1(heavyU1);
00063 heavyTerm += Gamma15(heavyU2)-Gamma15(heavyU1);
00064 heavyTerm /= 3.*std::sqrt(tm*EF);
00065 }
00066
00067 result = 0.5*(lightTerm+heavyTerm);
00068
00069 return result;
00070 }
00071
00072 G4double G4NeutronHPMadlandNixSpectrum::Sample(G4double anEnergy)
00073 {
00074 G4double tm = theMaxTemp.GetY(anEnergy);
00075 G4double last=0, buff, current = 100*MeV;
00076 G4double precision = 0.001;
00077 G4double newValue = 0., oldValue=0.;
00078 G4double random = G4UniformRand();
00079
00080 do
00081 {
00082 oldValue = newValue;
00083 newValue = FissionIntegral(tm, current);
00084 if(newValue < random)
00085 {
00086 buff = current;
00087 current+=std::abs(current-last)/2.;
00088 last = buff;
00089 if(current>190*MeV) throw G4HadronicException(__FILE__, __LINE__, "Madland-Nix Spectrum has not converged in sampling");
00090 }
00091 else
00092 {
00093 buff = current;
00094 current-=std::abs(current-last)/2.;
00095 last = buff;
00096 }
00097 }
00098 while (std::abs(oldValue-newValue)>precision*newValue);
00099 return current;
00100 }
00101
00102 G4double G4NeutronHPMadlandNixSpectrum::
00103 GIntegral(G4double tm, G4double anEnergy, G4double aMean)
00104 {
00105 if(aMean<1*eV) return 0;
00106 G4double b = anEnergy/eV;
00107 G4double sb = std::sqrt(b);
00108 G4double EF = aMean/eV;
00109
00110 G4double alpha = std::sqrt(tm);
00111 G4double beta = std::sqrt(EF);
00112 G4double A = EF/tm;
00113 G4double B = (sb+beta)*(sb+beta)/tm;
00114 G4double Ap = A;
00115 G4double Bp = (sb-beta)*(sb-beta)/tm;
00116
00117 G4double result;
00118 G4double alpha2 = alpha*alpha;
00119 G4double alphabeta = alpha*beta;
00120 if(b<EF)
00121 {
00122 result =
00123 (
00124 (0.4*alpha2*std::pow(B,2.5) - 0.5*alphabeta*B*B)*E1(B) -
00125 (0.4*alpha2*std::pow(A,2.5) - 0.5*alphabeta*A*A)*E1(A)
00126 )
00127 -
00128 (
00129 (0.4*alpha2*std::pow(Bp,2.5) + 0.5*alphabeta*Bp*Bp)*E1(Bp) -
00130 (0.4*alpha2*std::pow(Ap,2.5) + 0.5*alphabeta*Ap*Ap)*E1(Ap)
00131 )
00132 +
00133 (
00134 (alpha2*B-2*alphabeta*std::sqrt(B))*Gamma15(B) -
00135 (alpha2*A-2*alphabeta*std::sqrt(A))*Gamma15(A)
00136 )
00137 -
00138 (
00139 (alpha2*Bp-2*alphabeta*std::sqrt(Bp))*Gamma15(Bp) -
00140 (alpha2*Ap-2*alphabeta*std::sqrt(Ap))*Gamma15(Ap)
00141 )
00142 - 0.6*alpha2*(Gamma25(B) - Gamma25(A) - Gamma25(Bp) + Gamma25(Ap))
00143 - 1.5*alphabeta*(std::exp(-B)*(1+B) - std::exp(-A)*(1+A) + std::exp(-Bp)*(1+Bp) + std::exp(-Ap)*(1+Ap)) ;
00144 }
00145 else
00146 {
00147 result =
00148 (
00149 (0.4*alpha2*std::pow(B,2.5) - 0.5*alphabeta*B*B)*E1(B) -
00150 (0.4*alpha2*std::pow(A,2.5) - 0.5*alphabeta*A*A)*E1(A)
00151 );
00152 result -=
00153 (
00154 (0.4*alpha2*std::pow(Bp,2.5) + 0.5*alphabeta*Bp*Bp)*E1(Bp) -
00155 (0.4*alpha2*std::pow(Ap,2.5) + 0.5*alphabeta*Ap*Ap)*E1(Ap)
00156 );
00157 result +=
00158 (
00159 (alpha2*B-2*alphabeta*std::sqrt(B))*Gamma15(B) -
00160 (alpha2*A-2*alphabeta*std::sqrt(A))*Gamma15(A)
00161 );
00162 result -=
00163 (
00164 (alpha2*Bp+2*alphabeta*std::sqrt(Bp))*Gamma15(Bp) -
00165 (alpha2*Ap+2*alphabeta*std::sqrt(Ap))*Gamma15(Ap)
00166 );
00167 result -= 0.6*alpha2*(Gamma25(B) - Gamma25(A) - Gamma25(Bp) + Gamma25(Ap));
00168 result -= 1.5*alphabeta*(std::exp(-B)*(1+B) - std::exp(-A)*(1+A) + std::exp(-Bp)*(1+Bp) + std::exp(-Ap)*(1+Ap) - 2.) ;
00169 }
00170 result = result / (3.*std::sqrt(tm*EF));
00171 return result;
00172 }