Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NeutronHPMadlandNixSpectrum.hh
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 //
27 // $Id$
28 //
29 #ifndef G4NeutronHPMadlandNixSpectrum_h
30 #define G4NeutronHPMadlandNixSpectrum_h 1
31 
32 #include <fstream>
33 #include <cmath>
35 
36 #include "globals.hh"
37 #include "G4ios.hh"
38 #include "Randomize.hh"
39 #include "G4NeutronHPVector.hh"
40 #include "G4VNeutronHPEDis.hh"
41 
42 // #include <nag.h> @
43 // #include <nags.h> @
44 
45 
46 // we will need a List of these .... one per term.
47 
49 {
50  public:
52  {
53  expm1 = std::exp(-1.);
54  }
56  {
57  }
58 
59  inline void Init(std::ifstream & aDataFile)
60  {
61  theFractionalProb.Init(aDataFile);
62  aDataFile>> theAvarageKineticPerNucleonForLightFragments;
63  theAvarageKineticPerNucleonForLightFragments*=CLHEP::eV;
64  aDataFile>> theAvarageKineticPerNucleonForHeavyFragments;
65  theAvarageKineticPerNucleonForHeavyFragments*=CLHEP::eV;
66  theMaxTemp.Init(aDataFile);
67  }
68 
70  {
71  return theFractionalProb.GetY(anEnergy);
72  }
73 
74  G4double Sample(G4double anEnergy);
75 
76  private:
77 
78  G4double Madland(G4double aSecEnergy, G4double tm);
79 
80  inline G4double FissionIntegral(G4double tm, G4double anEnergy)
81  {
82  return 0.5*( GIntegral(tm, anEnergy, theAvarageKineticPerNucleonForLightFragments)
83  +GIntegral(tm, anEnergy, theAvarageKineticPerNucleonForHeavyFragments) );
84  }
85 
86  G4double GIntegral(G4double tm, G4double anEnergy, G4double aMean);
87 
88  inline G4double Gamma05(G4double aValue)
89  {
90  G4double result;
91  // gamma(1.2,x*X) = std::sqrt(CLHEP::pi)*Erf(x)
92  G4double x = std::sqrt(aValue);
93  G4double t = 1./(1+0.47047*x);
94  result = 1- (0.3480242*t - 0.0958798*t*t + 0.7478556*t*t*t)*std::exp(-aValue); // @ check
95  result *= std::sqrt(CLHEP::pi);
96  return result;
97  }
98 
99  inline G4double Gamma15(G4double aValue)
100  {
101  G4double result;
102  // gamma(a+1, x) = a*gamma(a,x)-x**a*std::exp(-x)
103  result = 0.5*Gamma05(aValue) - std::sqrt(aValue)*std::exp(-aValue); // @ check
104  return result;
105  }
106 
107  inline G4double Gamma25(G4double aValue)
108  {
109  G4double result;
110  result = 1.5*Gamma15(aValue) - std::pow(aValue,1.5)*std::exp(aValue); // @ check
111  return result;
112  }
113 
114  inline G4double E1(G4double aValue)
115  {
116  // good only for rather low aValue @@@ replace by the corresponding NAG function for the
117  // exponential integral. (<5 seems ok.
118  G4double gamma = 0.577216;
119  G4double precision = 0.000001;
120  G4double result =-gamma - std::log(aValue);
121  G4double term = -aValue;
122  //110527TKDB Unnessary codes, Detected by gcc4.6 compiler
123  //G4double last;
124  G4int count = 1;
125  result -= term;
126  for(;;)
127  {
128  count++;
129  //110527TKDB Unnessary codes, Detected by gcc4.6 compiler
130  //last = result;
131  term = -term*aValue*(count-1)/(count*count);
132  result -=term;
133  if(std::fabs(term)/std::fabs(result)<precision) break;
134  }
135 // NagError *fail; @
136 // result = nag_exp_integral(aValue, fail); @
137  return result;
138  }
139 
140  private:
141 
142  G4double expm1;
143 
144  private:
145 
146  G4NeutronHPVector theFractionalProb;
147 
148  G4double theAvarageKineticPerNucleonForLightFragments;
149  G4double theAvarageKineticPerNucleonForHeavyFragments;
150 
151  G4NeutronHPVector theMaxTemp;
152 
153 };
154 
155 #endif