Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPMadlandNixSpectrum.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 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
28 //
29 #ifndef G4ParticleHPMadlandNixSpectrum_h
30 #define G4ParticleHPMadlandNixSpectrum_h 1
31 
32 #include <fstream>
33 #include <cmath>
35 
36 #include "globals.hh"
37 #include "G4ios.hh"
38 #include "Randomize.hh"
39 #include "G4Exp.hh"
40 #include "G4Log.hh"
41 #include "G4Pow.hh"
42 #include "G4ParticleHPVector.hh"
43 #include "G4VParticleHPEDis.hh"
44 
45 // #include <nag.h> @
46 // #include <nags.h> @
47 
48 
49 // we will need a List of these .... one per term.
50 
52 {
53  public:
55  {
56  expm1 = G4Exp(-1.);
57  theAvarageKineticPerNucleonForLightFragments = 0.0;
58  theAvarageKineticPerNucleonForHeavyFragments = 0.0;
59  }
61  {
62  }
63 
64  inline void Init(std::istream & aDataFile)
65  {
66  theFractionalProb.Init(aDataFile);
67  aDataFile>> theAvarageKineticPerNucleonForLightFragments;
68  theAvarageKineticPerNucleonForLightFragments*=CLHEP::eV;
69  aDataFile>> theAvarageKineticPerNucleonForHeavyFragments;
70  theAvarageKineticPerNucleonForHeavyFragments*=CLHEP::eV;
71  theMaxTemp.Init(aDataFile);
72  }
73 
75  {
76  return theFractionalProb.GetY(anEnergy);
77  }
78 
79  G4double Sample(G4double anEnergy);
80 
81  private:
82 
83  G4double Madland(G4double aSecEnergy, G4double tm);
84 
85  inline G4double FissionIntegral(G4double tm, G4double anEnergy)
86  {
87  return 0.5*( GIntegral(tm, anEnergy, theAvarageKineticPerNucleonForLightFragments)
88  +GIntegral(tm, anEnergy, theAvarageKineticPerNucleonForHeavyFragments) );
89  }
90 
91  G4double GIntegral(G4double tm, G4double anEnergy, G4double aMean);
92 
93  inline G4double Gamma05(G4double aValue)
94  {
96  // gamma(1.2,x*X) = std::sqrt(CLHEP::pi)*Erf(x)
97  G4double x = std::sqrt(aValue);
98  G4double t = 1./(1+0.47047*x);
99  result = 1- (0.3480242*t - 0.0958798*t*t + 0.7478556*t*t*t)*G4Exp(-aValue); // @ check
100  result *= std::sqrt(CLHEP::pi);
101  return result;
102  }
103 
104  inline G4double Gamma15(G4double aValue)
105  {
107  // gamma(a+1, x) = a*gamma(a,x)-x**a*std::exp(-x)
108  result = 0.5*Gamma05(aValue) - std::sqrt(aValue)*G4Exp(-aValue); // @ check
109  return result;
110  }
111 
112  inline G4double Gamma25(G4double aValue)
113  {
115  result = 1.5*Gamma15(aValue) - G4Pow::GetInstance()->powA(aValue,1.5)*G4Exp(aValue); // @ check
116  return result;
117  }
118 
119  inline G4double E1(G4double aValue)
120  {
121  // good only for rather low aValue @@@ replace by the corresponding NAG function for the
122  // exponential integral. (<5 seems ok.
123  G4double gamma = 0.577216;
124  G4double precision = 0.000001;
125  G4double result =-gamma - G4Log(aValue);
126  G4double term = -aValue;
127  //110527TKDB Unnessary codes, Detected by gcc4.6 compiler
128  //G4double last;
129  G4int count = 1;
130  result -= term;
131  for(;;)
132  {
133  count++;
134  //110527TKDB Unnessary codes, Detected by gcc4.6 compiler
135  //last = result;
136  term = -term*aValue*(count-1)/(count*count);
137  result -=term;
138  if(std::fabs(term)/std::fabs(result)<precision) break;
139  }
140 // NagError *fail; @
141 // result = nag_exp_integral(aValue, fail); @
142  return result;
143  }
144 
145  private:
146 
147  G4double expm1;
148 
149  private:
150 
151  G4ParticleHPVector theFractionalProb;
152 
153  G4double theAvarageKineticPerNucleonForLightFragments;
154  G4double theAvarageKineticPerNucleonForHeavyFragments;
155 
156  G4ParticleHPVector theMaxTemp;
157 
158 };
159 
160 #endif
G4double G4ParticleHPJENDLHEData::G4double result
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
G4double GetFractionalProbability(G4double anEnergy)
void Init(std::istream &aDataFile, G4int total, G4double ux=1., G4double uy=1.)
G4double GetY(G4double x)
static constexpr double eV
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54