| Geant4 Cross Reference |
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 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
30 //
31 #include "G4ParticleHPMadlandNixSpectrum.hh"
32
33 #include "G4SystemOfUnits.hh"
34
35 G4double G4ParticleHPMadlandNixSpectrum::Madland(G4double aSecEnergy, G4double tm)
36 {
37 G4Pow* Pow = G4Pow::GetInstance();
38 G4double result;
39 G4double energy = aSecEnergy / eV;
40 G4double EF;
41
42 EF = theAvarageKineticPerNucleonForLightFragments / eV;
43 G4double lightU1 = std::sqrt(energy) - std::sqrt(EF);
44 lightU1 *= lightU1 / tm;
45 G4double lightU2 = std::sqrt(energy) + std::sqrt(EF);
46 lightU2 *= lightU2 / tm;
47 G4double lightTerm = 0;
48 if (theAvarageKineticPerNucleonForLightFragments > 1 * eV) {
49 lightTerm = Pow->powA(lightU2, 1.5) * E1(lightU2);
50 lightTerm -= Pow->powA(lightU1, 1.5) * E1(lightU1);
51 lightTerm += Gamma15(lightU2) - Gamma15(lightU1);
52 lightTerm /= 3. * std::sqrt(tm * EF);
53 }
54
55 EF = theAvarageKineticPerNucleonForHeavyFragments / eV;
56 G4double heavyU1 = std::sqrt(energy) - std::sqrt(EF);
57 heavyU1 *= heavyU1 / tm;
58 G4double heavyU2 = std::sqrt(energy) + std::sqrt(EF);
59 heavyU2 *= heavyU2 / tm;
60 G4double heavyTerm = 0;
61 if (theAvarageKineticPerNucleonForHeavyFragments > 1 * eV) {
62 heavyTerm = Pow->powA(heavyU2, 1.5) * E1(heavyU2);
63 heavyTerm -= Pow->powA(heavyU1, 1.5) * E1(heavyU1);
64 heavyTerm += Gamma15(heavyU2) - Gamma15(heavyU1);
65 heavyTerm /= 3. * std::sqrt(tm * EF);
66 }
67
68 result = 0.5 * (lightTerm + heavyTerm);
69
70 return result;
71 }
72
73 G4double G4ParticleHPMadlandNixSpectrum::Sample(G4double anEnergy)
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 threshold value at " << __LINE__ << "th line of "
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__, __LINE__,
98 "Madland-Nix Spectrum has not converged in sampling");
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 checking, 11.05.2015, T. Koi
107 return current;
108 }
109
110 G4double G4ParticleHPMadlandNixSpectrum::GIntegral(G4double tm, G4double anEnergy, G4double aMean)
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) - 0.5 * alphabeta * B * B) * E1(B)
130 - (0.4 * alpha2 * Pow->powA(A, 2.5) - 0.5 * alphabeta * A * A) * E1(A))
131 - ((0.4 * alpha2 * Pow->powA(Bp, 2.5) + 0.5 * alphabeta * Bp * Bp) * E1(Bp)
132 - (0.4 * alpha2 * Pow->powA(Ap, 2.5) + 0.5 * alphabeta * Ap * Ap) * E1(Ap))
133 + ((alpha2 * B - 2 * alphabeta * std::sqrt(B)) * Gamma15(B)
134 - (alpha2 * A - 2 * alphabeta * std::sqrt(A)) * Gamma15(A))
135 - ((alpha2 * Bp - 2 * alphabeta * std::sqrt(Bp)) * Gamma15(Bp)
136 - (alpha2 * Ap - 2 * alphabeta * std::sqrt(Ap)) * Gamma15(Ap))
137 - 0.6 * alpha2 * (Gamma25(B) - Gamma25(A) - Gamma25(Bp) + Gamma25(Ap))
138 - 1.5 * alphabeta
139 * (G4Exp(-B) * (1 + B) - G4Exp(-A) * (1 + A) + G4Exp(-Bp) * (1 + Bp)
140 + G4Exp(-Ap) * (1 + Ap));
141 }
142 else {
143 result = ((0.4 * alpha2 * Pow->powA(B, 2.5) - 0.5 * alphabeta * B * B) * E1(B)
144 - (0.4 * alpha2 * Pow->powA(A, 2.5) - 0.5 * alphabeta * A * A) * E1(A));
145 result -= ((0.4 * alpha2 * Pow->powA(Bp, 2.5) + 0.5 * alphabeta * Bp * Bp) * E1(Bp)
146 - (0.4 * alpha2 * Pow->powA(Ap, 2.5) + 0.5 * alphabeta * Ap * Ap) * E1(Ap));
147 result += ((alpha2 * B - 2 * alphabeta * std::sqrt(B)) * Gamma15(B)
148 - (alpha2 * A - 2 * alphabeta * std::sqrt(A)) * Gamma15(A));
149 result -= ((alpha2 * Bp + 2 * alphabeta * std::sqrt(Bp)) * Gamma15(Bp)
150 - (alpha2 * Ap + 2 * alphabeta * std::sqrt(Ap)) * Gamma15(Ap));
151 result -= 0.6 * alpha2 * (Gamma25(B) - Gamma25(A) - Gamma25(Bp) + Gamma25(Ap));
152 result -= 1.5 * alphabeta
153 * (G4Exp(-B) * (1 + B) - G4Exp(-A) * (1 + A) + G4Exp(-Bp) * (1 + Bp)
154 + G4Exp(-Ap) * (1 + Ap) - 2.);
155 }
156 result = result / (3. * std::sqrt(tm * EF));
157 return result;
158 }
159