41 G4double lightU1 = std::sqrt(energy)-std::sqrt(EF);
42 lightU1 *= lightU1/
tm;
43 G4double lightU2 = std::sqrt(energy)+std::sqrt(EF);
44 lightU2 *= lightU2/
tm;
48 lightTerm = std::pow(lightU2, 1.5)*
E1(lightU2);
49 lightTerm -= std::pow(lightU1, 1.5)*
E1(lightU1);
51 lightTerm /= 3.*std::sqrt(tm*EF);
55 G4double heavyU1 = std::sqrt(energy)-std::sqrt(EF);
56 heavyU1 *= heavyU1/
tm;
57 G4double heavyU2 = std::sqrt(energy)+std::sqrt(EF);
58 heavyU2 *= heavyU2/
tm;
62 heavyTerm = std::pow(heavyU2, 1.5)*
E1(heavyU2);
63 heavyTerm -= std::pow(heavyU1, 1.5)*
E1(heavyU1);
65 heavyTerm /= 3.*std::sqrt(tm*EF);
68 result = 0.5*(lightTerm+heavyTerm);
88 current+=std::abs(current-last)/2.;
90 if(current>190*
MeV)
throw G4HadronicException(__FILE__, __LINE__,
"Madland-Nix Spectrum has not converged in sampling");
95 current-=std::abs(current-last)/2.;
99 while (std::abs(oldValue-newValue)>precision*newValue);
106 if(aMean<1*
eV)
return 0;
125 (0.4*alpha2*std::pow(B,2.5) - 0.5*alphabeta*B*B)*
E1(B) -
126 (0.4*alpha2*std::pow(A,2.5) - 0.5*alphabeta*A*
A)*
E1(A)
130 (0.4*alpha2*std::pow(Bp,2.5) + 0.5*alphabeta*Bp*Bp)*
E1(Bp) -
131 (0.4*alpha2*std::pow(Ap,2.5) + 0.5*alphabeta*Ap*Ap)*
E1(Ap)
135 (alpha2*B-2*alphabeta*std::sqrt(B))*
Gamma15(B) -
136 (alpha2*A-2*alphabeta*std::sqrt(A))*
Gamma15(A)
140 (alpha2*Bp-2*alphabeta*std::sqrt(Bp))*
Gamma15(Bp) -
141 (alpha2*Ap-2*alphabeta*std::sqrt(Ap))*
Gamma15(Ap)
144 - 1.5*alphabeta*(std::exp(-B)*(1+B) - std::exp(-A)*(1+
A) + std::exp(-Bp)*(1+Bp) + std::exp(-Ap)*(1+Ap)) ;
150 (0.4*alpha2*std::pow(B,2.5) - 0.5*alphabeta*B*B)*
E1(B) -
151 (0.4*alpha2*std::pow(A,2.5) - 0.5*alphabeta*A*
A)*
E1(A)
155 (0.4*alpha2*std::pow(Bp,2.5) + 0.5*alphabeta*Bp*Bp)*
E1(Bp) -
156 (0.4*alpha2*std::pow(Ap,2.5) + 0.5*alphabeta*Ap*Ap)*
E1(Ap)
160 (alpha2*B-2*alphabeta*std::sqrt(B))*
Gamma15(B) -
161 (alpha2*A-2*alphabeta*std::sqrt(A))*
Gamma15(A)
165 (alpha2*Bp+2*alphabeta*std::sqrt(Bp))*
Gamma15(Bp) -
166 (alpha2*Ap+2*alphabeta*std::sqrt(Ap))*
Gamma15(Ap)
169 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.) ;
171 result = result / (3.*std::sqrt(tm*EF));
G4double theAvarageKineticPerNucleonForHeavyFragments
G4double Sample(G4double anEnergy)
G4double E1(G4double aValue)
G4double theAvarageKineticPerNucleonForLightFragments
G4double Gamma25(G4double aValue)
G4double GetY(G4double x)
static const G4double A[nN]
G4double Madland(G4double aSecEnergy, G4double tm)
G4double energy(const ThreeVector &p, const G4double m)
G4double Gamma15(G4double aValue)
G4ParticleHPVector theMaxTemp
G4double GIntegral(G4double tm, G4double anEnergy, G4double aMean)
G4double FissionIntegral(G4double tm, G4double anEnergy)
static const G4double alpha