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