source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPMadlandNixSpectrum.cc@ 1201

Last change on this file since 1201 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 6.1 KB
Line 
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 }
Note: See TracBrowser for help on using the repository browser.