Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPMadlandNixSpectrum.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 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
30 //
32 #include "G4SystemOfUnits.hh"
33 
34  G4double G4ParticleHPMadlandNixSpectrum::Madland(G4double aSecEnergy, G4double tm)
35  {
38  G4double energy = aSecEnergy/eV;
39  G4double EF;
40 
41  EF = theAvarageKineticPerNucleonForLightFragments/eV;
42  G4double lightU1 = std::sqrt(energy)-std::sqrt(EF);
43  lightU1 *= lightU1/tm;
44  G4double lightU2 = std::sqrt(energy)+std::sqrt(EF);
45  lightU2 *= lightU2/tm;
46  G4double lightTerm=0;
47  if(theAvarageKineticPerNucleonForLightFragments>1*eV)
48  {
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  {
63  heavyTerm = Pow->powA(heavyU2, 1.5)*E1(heavyU2);
64  heavyTerm -= Pow->powA(heavyU1, 1.5)*E1(heavyU1);
65  heavyTerm += Gamma15(heavyU2)-Gamma15(heavyU1);
66  heavyTerm /= 3.*std::sqrt(tm*EF);
67  }
68 
69  result = 0.5*(lightTerm+heavyTerm);
70 
71  return result;
72  }
73 
75  {
76  G4double tm = theMaxTemp.GetY(anEnergy);
77  G4double last=0, buff, current = 100*MeV;
78  G4double precision = 0.001;
79  G4double newValue = 0., oldValue=0.;
80  G4double random = G4UniformRand();
81 
82  G4int icounter=0;
83  G4int icounter_max=1024;
84  do
85  {
86  icounter++;
87  if ( icounter > icounter_max ) {
88  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
89  break;
90  }
91  oldValue = newValue;
92  newValue = FissionIntegral(tm, current);
93  if(newValue < random)
94  {
95  buff = current;
96  current+=std::abs(current-last)/2.;
97  last = buff;
98  if(current>190*MeV) throw G4HadronicException(__FILE__, __LINE__, "Madland-Nix Spectrum has not converged in sampling");
99  }
100  else
101  {
102  buff = current;
103  current-=std::abs(current-last)/2.;
104  last = buff;
105  }
106  }
107  while (std::abs(oldValue-newValue)>precision*newValue); // Loop checking, 11.05.2015, T. Koi
108  return current;
109  }
110 
111  G4double G4ParticleHPMadlandNixSpectrum::
112  GIntegral(G4double tm, G4double anEnergy, G4double aMean)
113  {
114  G4Pow* Pow=G4Pow::GetInstance();
115  if(aMean<1*eV) return 0;
116  G4double b = anEnergy/eV;
117  G4double sb = std::sqrt(b);
118  G4double EF = aMean/eV;
119 
120  G4double alpha = std::sqrt(tm);
121  G4double beta = std::sqrt(EF);
122  G4double A = EF/tm;
123  G4double B = (sb+beta)*(sb+beta)/tm;
124  G4double Ap = A;
125  G4double Bp = (sb-beta)*(sb-beta)/tm;
126 
128  G4double alpha2 = alpha*alpha;
129  G4double alphabeta = alpha*beta;
130  if(b<EF)
131  {
132  result =
133  (
134  (0.4*alpha2*Pow->powA(B,2.5) - 0.5*alphabeta*B*B)*E1(B) -
135  (0.4*alpha2*Pow->powA(A,2.5) - 0.5*alphabeta*A*A)*E1(A)
136  )
137  -
138  (
139  (0.4*alpha2*Pow->powA(Bp,2.5) + 0.5*alphabeta*Bp*Bp)*E1(Bp) -
140  (0.4*alpha2*Pow->powA(Ap,2.5) + 0.5*alphabeta*Ap*Ap)*E1(Ap)
141  )
142  +
143  (
144  (alpha2*B-2*alphabeta*std::sqrt(B))*Gamma15(B) -
145  (alpha2*A-2*alphabeta*std::sqrt(A))*Gamma15(A)
146  )
147  -
148  (
149  (alpha2*Bp-2*alphabeta*std::sqrt(Bp))*Gamma15(Bp) -
150  (alpha2*Ap-2*alphabeta*std::sqrt(Ap))*Gamma15(Ap)
151  )
152  - 0.6*alpha2*(Gamma25(B) - Gamma25(A) - Gamma25(Bp) + Gamma25(Ap))
153  - 1.5*alphabeta*(G4Exp(-B)*(1+B) - G4Exp(-A)*(1+A) + G4Exp(-Bp)*(1+Bp) + G4Exp(-Ap)*(1+Ap)) ;
154  }
155  else
156  {
157  result =
158  (
159  (0.4*alpha2*Pow->powA(B,2.5) - 0.5*alphabeta*B*B)*E1(B) -
160  (0.4*alpha2*Pow->powA(A,2.5) - 0.5*alphabeta*A*A)*E1(A)
161  );
162  result -=
163  (
164  (0.4*alpha2*Pow->powA(Bp,2.5) + 0.5*alphabeta*Bp*Bp)*E1(Bp) -
165  (0.4*alpha2*Pow->powA(Ap,2.5) + 0.5*alphabeta*Ap*Ap)*E1(Ap)
166  );
167  result +=
168  (
169  (alpha2*B-2*alphabeta*std::sqrt(B))*Gamma15(B) -
170  (alpha2*A-2*alphabeta*std::sqrt(A))*Gamma15(A)
171  );
172  result -=
173  (
174  (alpha2*Bp+2*alphabeta*std::sqrt(Bp))*Gamma15(Bp) -
175  (alpha2*Ap+2*alphabeta*std::sqrt(Ap))*Gamma15(Ap)
176  );
177  result -= 0.6*alpha2*(Gamma25(B) - Gamma25(A) - Gamma25(Bp) + Gamma25(Ap));
178  result -= 1.5*alphabeta*(G4Exp(-B)*(1+B) - G4Exp(-A)*(1+A) + G4Exp(-Bp)*(1+Bp) + G4Exp(-Ap)*(1+Ap) - 2.) ;
179  }
180  result = result / (3.*std::sqrt(tm*EF));
181  return result;
182  }
G4double G4ParticleHPJENDLHEData::G4double result
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
Definition: G4Pow.hh:56
double B(double temperature)
int G4int
Definition: G4Types.hh:78
tuple b
Definition: test.py:12
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
static constexpr double eV
Definition: G4SIunits.hh:215
G4double GetY(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double energy(const ThreeVector &p, const G4double m)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
double G4double
Definition: G4Types.hh:76
static const G4double alpha