| Geant4 Cross Reference |
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 "G4Exp.hh"
33 #include "G4Log.hh"
34 #include "G4ParticleHPVector.hh"
35 #include "G4Pow.hh"
36 #include "G4VParticleHPEDis.hh"
37 #include "G4ios.hh"
38 #include "Randomize.hh"
39 #include "globals.hh"
40
41 #include <CLHEP/Units/PhysicalConstants.h>
42
43 #include <cmath>
44 #include <fstream>
45
46 // #include <nag.h> @
47 // #include <nags.h> @
48
49 // we will need a List of these .... one per term.
50
51 class G4ParticleHPMadlandNixSpectrum : public G4VParticleHPEDis
52 {
53 public:
54 G4ParticleHPMadlandNixSpectrum()
55 {
56 expm1 = G4Exp(-1.);
57 theAvarageKineticPerNucleonForLightFragments = 0.0;
58 theAvarageKineticPerNucleonForHeavyFragments = 0.0;
59 }
60 ~G4ParticleHPMadlandNixSpectrum() override = default;
61
62 inline void Init(std::istream& aDataFile) override
63 {
64 theFractionalProb.Init(aDataFile);
65 aDataFile >> theAvarageKineticPerNucleonForLightFragments;
66 theAvarageKineticPerNucleonForLightFragments *= CLHEP::eV;
67 aDataFile >> theAvarageKineticPerNucleonForHeavyFragments;
68 theAvarageKineticPerNucleonForHeavyFragments *= CLHEP::eV;
69 theMaxTemp.Init(aDataFile);
70 }
71
72 inline G4double GetFractionalProbability(G4double anEnergy) override
73 {
74 return theFractionalProb.GetY(anEnergy);
75 }
76
77 G4double Sample(G4double anEnergy) override;
78
79 private:
80 G4double Madland(G4double aSecEnergy, G4double tm);
81
82 inline G4double FissionIntegral(G4double tm, G4double anEnergy)
83 {
84 return 0.5
85 * (GIntegral(tm, anEnergy, theAvarageKineticPerNucleonForLightFragments)
86 + GIntegral(tm, anEnergy, theAvarageKineticPerNucleonForHeavyFragments));
87 }
88
89 G4double GIntegral(G4double tm, G4double anEnergy, G4double aMean);
90
91 inline G4double Gamma05(G4double aValue)
92 {
93 G4double result;
94 // gamma(1.2,x*X) = std::sqrt(CLHEP::pi)*Erf(x)
95 G4double x = std::sqrt(aValue);
96 G4double t = 1. / (1 + 0.47047 * x);
97 result =
98 1
99 - (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 =
116 1.5 * Gamma15(aValue) - G4Pow::GetInstance()->powA(aValue, 1.5) * G4Exp(aValue); // @ check
117 return result;
118 }
119
120 inline G4double E1(G4double aValue)
121 {
122 // good only for rather low aValue @@@ replace by the corresponding NAG function for the
123 // exponential integral. (<5 seems ok.
124 G4double gamma = 0.577216;
125 G4double precision = 0.000001;
126 G4double result = -gamma - G4Log(aValue);
127 G4double term = -aValue;
128 // 110527TKDB Unnessary codes, Detected by gcc4.6 compiler
129 // G4double last;
130 G4int count = 1;
131 result -= term;
132 for (;;) {
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 G4double expm1;
147
148 private:
149 G4ParticleHPVector theFractionalProb;
150
151 G4double theAvarageKineticPerNucleonForLightFragments;
152 G4double theAvarageKineticPerNucleonForHeavyFragments;
153
154 G4ParticleHPVector theMaxTemp;
155 };
156
157 #endif
158