Geant4_10
G4NeutronHPMadlandNixSpectrum.cc
Go to the documentation of this file.
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 
31 #include "G4SystemOfUnits.hh"
32 
33  G4double G4NeutronHPMadlandNixSpectrum::Madland(G4double aSecEnergy, G4double tm)
34  {
36  G4double energy = aSecEnergy/eV;
37  G4double EF;
38 
39  EF = theAvarageKineticPerNucleonForLightFragments/eV;
40  G4double lightU1 = std::sqrt(energy)-std::sqrt(EF);
41  lightU1 *= lightU1/tm;
42  G4double lightU2 = std::sqrt(energy)+std::sqrt(EF);
43  lightU2 *= lightU2/tm;
44  G4double lightTerm=0;
45  if(theAvarageKineticPerNucleonForLightFragments>1*eV)
46  {
47  lightTerm = std::pow(lightU2, 1.5)*E1(lightU2);
48  lightTerm -= std::pow(lightU1, 1.5)*E1(lightU1);
49  lightTerm += Gamma15(lightU2)-Gamma15(lightU1);
50  lightTerm /= 3.*std::sqrt(tm*EF);
51  }
52 
53  EF = theAvarageKineticPerNucleonForHeavyFragments/eV;
54  G4double heavyU1 = std::sqrt(energy)-std::sqrt(EF);
55  heavyU1 *= heavyU1/tm;
56  G4double heavyU2 = std::sqrt(energy)+std::sqrt(EF);
57  heavyU2 *= heavyU2/tm;
58  G4double heavyTerm=0 ;
59  if(theAvarageKineticPerNucleonForHeavyFragments> 1*eV)
60  {
61  heavyTerm = std::pow(heavyU2, 1.5)*E1(heavyU2);
62  heavyTerm -= std::pow(heavyU1, 1.5)*E1(heavyU1);
63  heavyTerm += Gamma15(heavyU2)-Gamma15(heavyU1);
64  heavyTerm /= 3.*std::sqrt(tm*EF);
65  }
66 
67  result = 0.5*(lightTerm+heavyTerm);
68 
69  return result;
70  }
71 
73  {
74  G4double tm = theMaxTemp.GetY(anEnergy);
75  G4double last=0, buff, current = 100*MeV;
76  G4double precision = 0.001;
77  G4double newValue = 0., oldValue=0.;
78  G4double random = G4UniformRand();
79 
80  do
81  {
82  oldValue = newValue;
83  newValue = FissionIntegral(tm, current);
84  if(newValue < random)
85  {
86  buff = current;
87  current+=std::abs(current-last)/2.;
88  last = buff;
89  if(current>190*MeV) throw G4HadronicException(__FILE__, __LINE__, "Madland-Nix Spectrum has not converged in sampling");
90  }
91  else
92  {
93  buff = current;
94  current-=std::abs(current-last)/2.;
95  last = buff;
96  }
97  }
98  while (std::abs(oldValue-newValue)>precision*newValue);
99  return current;
100  }
101 
102  G4double G4NeutronHPMadlandNixSpectrum::
103  GIntegral(G4double tm, G4double anEnergy, G4double aMean)
104  {
105  if(aMean<1*eV) return 0;
106  G4double b = anEnergy/eV;
107  G4double sb = std::sqrt(b);
108  G4double EF = aMean/eV;
109 
110  G4double alpha = std::sqrt(tm);
111  G4double beta = std::sqrt(EF);
112  G4double A = EF/tm;
113  G4double B = (sb+beta)*(sb+beta)/tm;
114  G4double Ap = A;
115  G4double Bp = (sb-beta)*(sb-beta)/tm;
116 
118  G4double alpha2 = alpha*alpha;
119  G4double alphabeta = alpha*beta;
120  if(b<EF)
121  {
122  result =
123  (
124  (0.4*alpha2*std::pow(B,2.5) - 0.5*alphabeta*B*B)*E1(B) -
125  (0.4*alpha2*std::pow(A,2.5) - 0.5*alphabeta*A*A)*E1(A)
126  )
127  -
128  (
129  (0.4*alpha2*std::pow(Bp,2.5) + 0.5*alphabeta*Bp*Bp)*E1(Bp) -
130  (0.4*alpha2*std::pow(Ap,2.5) + 0.5*alphabeta*Ap*Ap)*E1(Ap)
131  )
132  +
133  (
134  (alpha2*B-2*alphabeta*std::sqrt(B))*Gamma15(B) -
135  (alpha2*A-2*alphabeta*std::sqrt(A))*Gamma15(A)
136  )
137  -
138  (
139  (alpha2*Bp-2*alphabeta*std::sqrt(Bp))*Gamma15(Bp) -
140  (alpha2*Ap-2*alphabeta*std::sqrt(Ap))*Gamma15(Ap)
141  )
142  - 0.6*alpha2*(Gamma25(B) - Gamma25(A) - Gamma25(Bp) + Gamma25(Ap))
143  - 1.5*alphabeta*(std::exp(-B)*(1+B) - std::exp(-A)*(1+A) + std::exp(-Bp)*(1+Bp) + std::exp(-Ap)*(1+Ap)) ;
144  }
145  else
146  {
147  result =
148  (
149  (0.4*alpha2*std::pow(B,2.5) - 0.5*alphabeta*B*B)*E1(B) -
150  (0.4*alpha2*std::pow(A,2.5) - 0.5*alphabeta*A*A)*E1(A)
151  );
152  result -=
153  (
154  (0.4*alpha2*std::pow(Bp,2.5) + 0.5*alphabeta*Bp*Bp)*E1(Bp) -
155  (0.4*alpha2*std::pow(Ap,2.5) + 0.5*alphabeta*Ap*Ap)*E1(Ap)
156  );
157  result +=
158  (
159  (alpha2*B-2*alphabeta*std::sqrt(B))*Gamma15(B) -
160  (alpha2*A-2*alphabeta*std::sqrt(A))*Gamma15(A)
161  );
162  result -=
163  (
164  (alpha2*Bp+2*alphabeta*std::sqrt(Bp))*Gamma15(Bp) -
165  (alpha2*Ap+2*alphabeta*std::sqrt(Ap))*Gamma15(Ap)
166  );
167  result -= 0.6*alpha2*(Gamma25(B) - Gamma25(A) - Gamma25(Bp) + Gamma25(Ap));
168  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.) ;
169  }
170  result = result / (3.*std::sqrt(tm*EF));
171  return result;
172  }
G4double GetY(G4double x)
G4double G4NeutronHPJENDLHEData::G4double result
tuple b
Definition: test.py:12
#define G4UniformRand()
Definition: Randomize.hh:87
double energy
Definition: plottest35.C:25
double G4double
Definition: G4Types.hh:76