| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans
29 // P. Arce, June-2014 Conversion neutron_hp to
30 //
31 #include "G4ParticleHPMadlandNixSpectrum.hh"
32
33 #include "G4SystemOfUnits.hh"
34
35 G4double G4ParticleHPMadlandNixSpectrum::Madla
36 {
37 G4Pow* Pow = G4Pow::GetInstance();
38 G4double result;
39 G4double energy = aSecEnergy / eV;
40 G4double EF;
41
42 EF = theAvarageKineticPerNucleonForLightFrag
43 G4double lightU1 = std::sqrt(energy) - std::
44 lightU1 *= lightU1 / tm;
45 G4double lightU2 = std::sqrt(energy) + std::
46 lightU2 *= lightU2 / tm;
47 G4double lightTerm = 0;
48 if (theAvarageKineticPerNucleonForLightFragm
49 lightTerm = Pow->powA(lightU2, 1.5) * E1(l
50 lightTerm -= Pow->powA(lightU1, 1.5) * E1(
51 lightTerm += Gamma15(lightU2) - Gamma15(li
52 lightTerm /= 3. * std::sqrt(tm * EF);
53 }
54
55 EF = theAvarageKineticPerNucleonForHeavyFrag
56 G4double heavyU1 = std::sqrt(energy) - std::
57 heavyU1 *= heavyU1 / tm;
58 G4double heavyU2 = std::sqrt(energy) + std::
59 heavyU2 *= heavyU2 / tm;
60 G4double heavyTerm = 0;
61 if (theAvarageKineticPerNucleonForHeavyFragm
62 heavyTerm = Pow->powA(heavyU2, 1.5) * E1(h
63 heavyTerm -= Pow->powA(heavyU1, 1.5) * E1(
64 heavyTerm += Gamma15(heavyU2) - Gamma15(he
65 heavyTerm /= 3. * std::sqrt(tm * EF);
66 }
67
68 result = 0.5 * (lightTerm + heavyTerm);
69
70 return result;
71 }
72
73 G4double G4ParticleHPMadlandNixSpectrum::Sampl
74 {
75 G4double tm = theMaxTemp.GetY(anEnergy);
76 G4double last = 0, buff, current = 100 * MeV
77 G4double precision = 0.001;
78 G4double newValue = 0., oldValue = 0.;
79 G4double random = G4UniformRand();
80
81 G4int icounter = 0;
82 G4int icounter_max = 1024;
83 do {
84 icounter++;
85 if (icounter > icounter_max) {
86 G4cout << "Loop-counter exceeded the thr
87 << __FILE__ << "." << G4endl;
88 break;
89 }
90 oldValue = newValue;
91 newValue = FissionIntegral(tm, current);
92 if (newValue < random) {
93 buff = current;
94 current += std::abs(current - last) / 2.
95 last = buff;
96 if (current > 190 * MeV)
97 throw G4HadronicException(__FILE__, __
98 "Madland-Nix
99 }
100 else {
101 buff = current;
102 current -= std::abs(current - last) / 2.
103 last = buff;
104 }
105 } while (std::abs(oldValue - newValue)
106 > precision * newValue); // Loop c
107 return current;
108 }
109
110 G4double G4ParticleHPMadlandNixSpectrum::GInte
111 {
112 G4Pow* Pow = G4Pow::GetInstance();
113 if (aMean < 1 * eV) return 0;
114 G4double b = anEnergy / eV;
115 G4double sb = std::sqrt(b);
116 G4double EF = aMean / eV;
117
118 G4double alpha = std::sqrt(tm);
119 G4double beta = std::sqrt(EF);
120 G4double A = EF / tm;
121 G4double B = (sb + beta) * (sb + beta) / tm;
122 G4double Ap = A;
123 G4double Bp = (sb - beta) * (sb - beta) / tm
124
125 G4double result;
126 G4double alpha2 = alpha * alpha;
127 G4double alphabeta = alpha * beta;
128 if (b < EF) {
129 result = ((0.4 * alpha2 * Pow->powA(B, 2.5
130 - (0.4 * alpha2 * Pow->powA(A, 2
131 - ((0.4 * alpha2 * Pow->powA(Bp,
132 - (0.4 * alpha2 * Pow->powA(Ap
133 + ((alpha2 * B - 2 * alphabeta *
134 - (alpha2 * A - 2 * alphabeta
135 - ((alpha2 * Bp - 2 * alphabeta *
136 - (alpha2 * Ap - 2 * alphabeta
137 - 0.6 * alpha2 * (Gamma25(B) - Ga
138 - 1.5 * alphabeta
139 * (G4Exp(-B) * (1 + B) - G4Ex
140 + G4Exp(-Ap) * (1 + Ap));
141 }
142 else {
143 result = ((0.4 * alpha2 * Pow->powA(B, 2.5
144 - (0.4 * alpha2 * Pow->powA(A, 2
145 result -= ((0.4 * alpha2 * Pow->powA(Bp, 2
146 - (0.4 * alpha2 * Pow->powA(Ap,
147 result += ((alpha2 * B - 2 * alphabeta * s
148 - (alpha2 * A - 2 * alphabeta *
149 result -= ((alpha2 * Bp + 2 * alphabeta *
150 - (alpha2 * Ap + 2 * alphabeta
151 result -= 0.6 * alpha2 * (Gamma25(B) - Gam
152 result -= 1.5 * alphabeta
153 * (G4Exp(-B) * (1 + B) - G4Exp(-
154 + G4Exp(-Ap) * (1 + Ap) - 2.)
155 }
156 result = result / (3. * std::sqrt(tm * EF));
157 return result;
158 }
159