ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPMadlandNixSpectrum.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPMadlandNixSpectrum.hh
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.);
59  }
61  {
62  }
63 
64  inline void Init(std::istream & aDataFile)
65  {
66  theFractionalProb.Init(aDataFile);
68  theAvarageKineticPerNucleonForLightFragments*=CLHEP::eV;
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 
86  {
87  return 0.5*( GIntegral(tm, anEnergy, theAvarageKineticPerNucleonForLightFragments)
89  }
90 
91  G4double GIntegral(G4double tm, G4double anEnergy, G4double aMean);
92 
93  inline G4double Gamma05(G4double aValue)
94  {
95  G4double result;
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  {
106  G4double result;
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  {
114  G4double result;
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 
148 
149  private:
150 
152 
155 
157 
158 };
159 
160 #endif