1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // neutron_hp -- source file |
---|
27 | // J.P. Wellisch, Nov-1996 |
---|
28 | // A prototype of the low energy neutron transport model. |
---|
29 | #include "G4NeutronHPMadlandNixSpectrum.hh" |
---|
30 | |
---|
31 | G4double G4NeutronHPMadlandNixSpectrum::Madland(G4double aSecEnergy, G4double tm) |
---|
32 | { |
---|
33 | G4double result; |
---|
34 | G4double energy = aSecEnergy/eV; |
---|
35 | G4double EF; |
---|
36 | |
---|
37 | EF = theAvarageKineticPerNucleonForLightFragments/eV; |
---|
38 | G4double lightU1 = std::sqrt(energy)-std::sqrt(EF); |
---|
39 | lightU1 *= lightU1/tm; |
---|
40 | G4double lightU2 = std::sqrt(energy)+std::sqrt(EF); |
---|
41 | lightU2 *= lightU2/tm; |
---|
42 | G4double lightTerm=0; |
---|
43 | if(theAvarageKineticPerNucleonForLightFragments>1*eV) |
---|
44 | { |
---|
45 | lightTerm = std::pow(lightU2, 1.5)*E1(lightU2); |
---|
46 | lightTerm -= std::pow(lightU1, 1.5)*E1(lightU1); |
---|
47 | lightTerm += Gamma15(lightU2)-Gamma15(lightU1); |
---|
48 | lightTerm /= 3.*std::sqrt(tm*EF); |
---|
49 | } |
---|
50 | |
---|
51 | EF = theAvarageKineticPerNucleonForHeavyFragments/eV; |
---|
52 | G4double heavyU1 = std::sqrt(energy)-std::sqrt(EF); |
---|
53 | heavyU1 *= heavyU1/tm; |
---|
54 | G4double heavyU2 = std::sqrt(energy)+std::sqrt(EF); |
---|
55 | heavyU2 *= heavyU2/tm; |
---|
56 | G4double heavyTerm=0 ; |
---|
57 | if(theAvarageKineticPerNucleonForHeavyFragments> 1*eV) |
---|
58 | { |
---|
59 | heavyTerm = std::pow(heavyU2, 1.5)*E1(heavyU2); |
---|
60 | heavyTerm -= std::pow(heavyU1, 1.5)*E1(heavyU1); |
---|
61 | heavyTerm += Gamma15(heavyU2)-Gamma15(heavyU1); |
---|
62 | heavyTerm /= 3.*std::sqrt(tm*EF); |
---|
63 | } |
---|
64 | |
---|
65 | result = 0.5*(lightTerm+heavyTerm); |
---|
66 | |
---|
67 | return result; |
---|
68 | } |
---|
69 | |
---|
70 | G4double G4NeutronHPMadlandNixSpectrum::Sample(G4double anEnergy) |
---|
71 | { |
---|
72 | G4double tm = theMaxTemp.GetY(anEnergy); |
---|
73 | G4double last=0, buff, current = 100*MeV; |
---|
74 | G4double precision = 0.001; |
---|
75 | G4double newValue = 0., oldValue=0.; |
---|
76 | G4double random = G4UniformRand(); |
---|
77 | |
---|
78 | do |
---|
79 | { |
---|
80 | oldValue = newValue; |
---|
81 | newValue = FissionIntegral(tm, current); |
---|
82 | if(newValue < random) |
---|
83 | { |
---|
84 | buff = current; |
---|
85 | current+=std::abs(current-last)/2.; |
---|
86 | last = buff; |
---|
87 | if(current>190*MeV) throw G4HadronicException(__FILE__, __LINE__, "Madland-Nix Spectrum has not converged in sampling"); |
---|
88 | } |
---|
89 | else |
---|
90 | { |
---|
91 | buff = current; |
---|
92 | current-=std::abs(current-last)/2.; |
---|
93 | last = buff; |
---|
94 | } |
---|
95 | } |
---|
96 | while (std::abs(oldValue-newValue)>precision*newValue); |
---|
97 | return current; |
---|
98 | } |
---|
99 | |
---|
100 | G4double G4NeutronHPMadlandNixSpectrum:: |
---|
101 | GIntegral(G4double tm, G4double anEnergy, G4double aMean) |
---|
102 | { |
---|
103 | if(aMean<1*eV) return 0; |
---|
104 | G4double b = anEnergy/eV; |
---|
105 | G4double sb = std::sqrt(b); |
---|
106 | G4double EF = aMean/eV; |
---|
107 | |
---|
108 | G4double alpha = std::sqrt(tm); |
---|
109 | G4double beta = std::sqrt(EF); |
---|
110 | G4double A = EF/tm; |
---|
111 | G4double B = (sb+beta)*(sb+beta)/tm; |
---|
112 | G4double Ap = A; |
---|
113 | G4double Bp = (sb-beta)*(sb-beta)/tm; |
---|
114 | |
---|
115 | G4double result; |
---|
116 | G4double alpha2 = alpha*alpha; |
---|
117 | G4double alphabeta = alpha*beta; |
---|
118 | if(b<EF) |
---|
119 | { |
---|
120 | result = |
---|
121 | ( |
---|
122 | (0.4*alpha2*std::pow(B,2.5) - 0.5*alphabeta*B*B)*E1(B) - |
---|
123 | (0.4*alpha2*std::pow(A,2.5) - 0.5*alphabeta*A*A)*E1(A) |
---|
124 | ) |
---|
125 | - |
---|
126 | ( |
---|
127 | (0.4*alpha2*std::pow(Bp,2.5) + 0.5*alphabeta*Bp*Bp)*E1(Bp) - |
---|
128 | (0.4*alpha2*std::pow(Ap,2.5) + 0.5*alphabeta*Ap*Ap)*E1(Ap) |
---|
129 | ) |
---|
130 | + |
---|
131 | ( |
---|
132 | (alpha2*B-2*alphabeta*std::sqrt(B))*Gamma15(B) - |
---|
133 | (alpha2*A-2*alphabeta*std::sqrt(A))*Gamma15(A) |
---|
134 | ) |
---|
135 | - |
---|
136 | ( |
---|
137 | (alpha2*Bp-2*alphabeta*std::sqrt(Bp))*Gamma15(Bp) - |
---|
138 | (alpha2*Ap-2*alphabeta*std::sqrt(Ap))*Gamma15(Ap) |
---|
139 | ) |
---|
140 | - 0.6*alpha2*(Gamma25(B) - Gamma25(A) - Gamma25(Bp) + Gamma25(Ap)) |
---|
141 | - 1.5*alphabeta*(std::exp(-B)*(1+B) - std::exp(-A)*(1+A) + std::exp(-Bp)*(1+Bp) + std::exp(-Ap)*(1+Ap)) ; |
---|
142 | } |
---|
143 | else |
---|
144 | { |
---|
145 | result = |
---|
146 | ( |
---|
147 | (0.4*alpha2*std::pow(B,2.5) - 0.5*alphabeta*B*B)*E1(B) - |
---|
148 | (0.4*alpha2*std::pow(A,2.5) - 0.5*alphabeta*A*A)*E1(A) |
---|
149 | ); |
---|
150 | result -= |
---|
151 | ( |
---|
152 | (0.4*alpha2*std::pow(Bp,2.5) + 0.5*alphabeta*Bp*Bp)*E1(Bp) - |
---|
153 | (0.4*alpha2*std::pow(Ap,2.5) + 0.5*alphabeta*Ap*Ap)*E1(Ap) |
---|
154 | ); |
---|
155 | result += |
---|
156 | ( |
---|
157 | (alpha2*B-2*alphabeta*std::sqrt(B))*Gamma15(B) - |
---|
158 | (alpha2*A-2*alphabeta*std::sqrt(A))*Gamma15(A) |
---|
159 | ); |
---|
160 | result -= |
---|
161 | ( |
---|
162 | (alpha2*Bp+2*alphabeta*std::sqrt(Bp))*Gamma15(Bp) - |
---|
163 | (alpha2*Ap+2*alphabeta*std::sqrt(Ap))*Gamma15(Ap) |
---|
164 | ); |
---|
165 | result -= 0.6*alpha2*(Gamma25(B) - Gamma25(A) - Gamma25(Bp) + Gamma25(Ap)); |
---|
166 | 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.) ; |
---|
167 | } |
---|
168 | result = result / (3.*std::sqrt(tm*EF)); |
---|
169 | return result; |
---|
170 | } |
---|