| 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 // 080809 Change interpolation scheme of "histogram", now using LinearLinear
28 // For multidimensional interpolations By T. Koi
29 //
30 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
31 //
32 #ifndef G4ParticleHPInterpolator_h
33 #define G4ParticleHPInterpolator_h 1
34
35 #include "G4Exp.hh"
36 #include "G4HadronicException.hh"
37 #include "G4InterpolationScheme.hh"
38 #include "G4Log.hh"
39 #include "G4ios.hh"
40 #include "Randomize.hh"
41 #include "globals.hh"
42
43 class G4ParticleHPInterpolator
44 {
45 public:
46 G4ParticleHPInterpolator() = default;
47 ~G4ParticleHPInterpolator() = default;
48
49 inline G4double Lin(G4double x, G4double x1, G4double x2, G4double y1, G4double y2)
50 {
51 G4double slope = 0, off = 0;
52 if (x2 - x1 == 0) return (y2 + y1) / 2.;
53 slope = (y2 - y1) / (x2 - x1);
54 off = y2 - x2 * slope;
55 G4double y = x * slope + off;
56 return y;
57 }
58
59 inline G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2,
60 G4double y1, G4double y2) const;
61 inline G4double Interpolate2(G4InterpolationScheme aScheme, G4double x, G4double x1,
62 G4double x2, G4double y1, G4double y2) const;
63
64 G4double GetBinIntegral(const G4InterpolationScheme& aScheme, const G4double x1,
65 const G4double x2, const G4double y1, const G4double y2);
66
67 G4double GetWeightedBinIntegral(const G4InterpolationScheme& aScheme, const G4double x1,
68 const G4double x2, const G4double y1, const G4double y2);
69
70 private:
71 inline G4double Histogram(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const;
72 inline G4double LinearLinear(G4double x, G4double x1, G4double x2, G4double y1,
73 G4double y2) const;
74 inline G4double LinearLogarithmic(G4double x, G4double x1, G4double x2, G4double y1,
75 G4double y2) const;
76 inline G4double LogarithmicLinear(G4double x, G4double x1, G4double x2, G4double y1,
77 G4double y2) const;
78 inline G4double LogarithmicLogarithmic(G4double x, G4double x1, G4double x2, G4double y1,
79 G4double y2) const;
80 inline G4double Random(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const;
81 };
82
83 inline G4double G4ParticleHPInterpolator::Interpolate(G4InterpolationScheme aScheme, G4double x,
84 G4double x1, G4double x2, G4double y1,
85 G4double y2) const
86 {
87 G4double result(0);
88 G4int theScheme = aScheme;
89 theScheme = theScheme % CSTART_;
90 switch (theScheme) {
91 case 1:
92 // 080809
93 // result = Histogram(x, x1, x2, y1, y2);
94 result = LinearLinear(x, x1, x2, y1, y2);
95 break;
96 case 2:
97 result = LinearLinear(x, x1, x2, y1, y2);
98 break;
99 case 3:
100 result = LinearLogarithmic(x, x1, x2, y1, y2);
101 break;
102 case 4:
103 result = LogarithmicLinear(x, x1, x2, y1, y2);
104 break;
105 case 5:
106 result = LogarithmicLogarithmic(x, x1, x2, y1, y2);
107 break;
108 case 6:
109 result = Random(x, x1, x2, y1, y2);
110 break;
111 default:
112 G4cout << "theScheme = " << theScheme << G4endl;
113 throw G4HadronicException(__FILE__, __LINE__,
114 "G4ParticleHPInterpolator::Carthesian Invalid InterpolationScheme");
115 break;
116 }
117 return result;
118 }
119
120 inline G4double G4ParticleHPInterpolator::Interpolate2(G4InterpolationScheme aScheme, G4double x,
121 G4double x1, G4double x2, G4double y1,
122 G4double y2) const
123 {
124 G4double result(0);
125 G4int theScheme = aScheme;
126 theScheme = theScheme % CSTART_;
127 switch (theScheme) {
128 case 1:
129 result = Histogram(x, x1, x2, y1, y2);
130 break;
131 case 2:
132 result = LinearLinear(x, x1, x2, y1, y2);
133 break;
134 case 3:
135 result = LinearLogarithmic(x, x1, x2, y1, y2);
136 break;
137 case 4:
138 result = LogarithmicLinear(x, x1, x2, y1, y2);
139 break;
140 case 5:
141 result = LogarithmicLogarithmic(x, x1, x2, y1, y2);
142 break;
143 case 6:
144 result = Random(x, x1, x2, y1, y2);
145 break;
146 default:
147 G4cout << "theScheme = " << theScheme << G4endl;
148 throw G4HadronicException(__FILE__, __LINE__,
149 "G4ParticleHPInterpolator::Carthesian Invalid InterpolationScheme");
150 break;
151 }
152 return result;
153 }
154
155 inline G4double G4ParticleHPInterpolator::Histogram(G4double, G4double, G4double, G4double y1,
156 G4double) const
157 {
158 G4double result;
159 result = y1;
160 return result;
161 }
162
163 inline G4double G4ParticleHPInterpolator::LinearLinear(G4double x, G4double x1, G4double x2,
164 G4double y1, G4double y2) const
165 {
166 G4double slope = 0, off = 0;
167 if (x2 - x1 == 0) return (y2 + y1) / 2.;
168 slope = (y2 - y1) / (x2 - x1);
169 off = y2 - x2 * slope;
170 G4double y = x * slope + off;
171 return y;
172 }
173
174 inline G4double G4ParticleHPInterpolator::LinearLogarithmic(G4double x, G4double x1, G4double x2,
175 G4double y1, G4double y2) const
176 {
177 G4double result;
178 if (x == 0)
179 result = y1 + y2 / 2.;
180 else if (x1 == 0)
181 result = y1;
182 else if (x2 == 0)
183 result = y2;
184 else
185 result = LinearLinear(G4Log(x), G4Log(x1), G4Log(x2), y1, y2);
186 return result;
187 }
188
189 inline G4double G4ParticleHPInterpolator::LogarithmicLinear(G4double x, G4double x1, G4double x2,
190 G4double y1, G4double y2) const
191 {
192 G4double result;
193 if (y1 == 0 || y2 == 0)
194 result = 0;
195 else {
196 result = LinearLinear(x, x1, x2, G4Log(y1), G4Log(y2));
197 result = G4Exp(result);
198 }
199 return result;
200 }
201
202 inline G4double G4ParticleHPInterpolator::LogarithmicLogarithmic(G4double x, G4double x1,
203 G4double x2, G4double y1,
204 G4double y2) const
205 {
206 if (x == 0) return y1 + y2 / 2.;
207 if (x1 == 0) return y1;
208 if (x2 == 0) return y2;
209 G4double result;
210 if (y1 == 0 || y2 == 0)
211 result = 0;
212 else {
213 result = LinearLinear(G4Log(x), G4Log(x1), G4Log(x2), G4Log(y1), G4Log(y2));
214 result = G4Exp(result);
215 }
216 return result;
217 }
218
219 inline G4double G4ParticleHPInterpolator::Random(G4double, G4double, G4double, G4double y1,
220 G4double y2) const
221 {
222 G4double result;
223 result = y1 + G4UniformRand() * (y2 - y1);
224 return result;
225 }
226
227 #endif
228