| 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 // G4PhysicsVector inline methods implementation
27 //
28 // Authors:
29 // - 02 Dec. 1995, G.Cosmo: Structure created based on object model
30 // - 03 Mar. 1996, K.Amako: Implemented the 1st version
31 // --------------------------------------------------------------------
32 inline G4double G4PhysicsVector::operator[](const std::size_t index) const
33 {
34 return dataVector[index];
35 }
36
37 // ---------------------------------------------------------------
38 inline G4double G4PhysicsVector::operator()(const std::size_t index) const
39 {
40 return dataVector[index];
41 }
42
43 // ---------------------------------------------------------------
44 inline G4double G4PhysicsVector::Energy(const std::size_t index) const
45 {
46 return binVector[index];
47 }
48
49 // ---------------------------------------------------------------
50 inline G4double
51 G4PhysicsVector::GetLowEdgeEnergy(const std::size_t index) const
52 {
53 return binVector[index];
54 }
55
56 // ---------------------------------------------------------------
57 inline G4double G4PhysicsVector::GetMinEnergy() const
58 {
59 return edgeMin;
60 }
61
62 // ---------------------------------------------------------------
63 inline G4double G4PhysicsVector::GetMaxEnergy() const
64 {
65 return edgeMax;
66 }
67
68 // ---------------------------------------------------------------
69 inline G4double G4PhysicsVector::GetMinValue() const
70 {
71 return (numberOfNodes > 0) ? dataVector[0] : 0.0;
72 }
73
74 // ---------------------------------------------------------------
75 inline G4double G4PhysicsVector::GetMaxValue() const
76 {
77 return (numberOfNodes > 0) ? dataVector[numberOfNodes - 1] : 0.0;
78 }
79
80 // ---------------------------------------------------------------
81 inline std::size_t G4PhysicsVector::GetVectorLength() const
82 {
83 return numberOfNodes;
84 }
85
86 // ---------------------------------------------------------------
87 inline void G4PhysicsVector::PutValue(std::size_t index, G4double theValue)
88 {
89 if(index >= numberOfNodes)
90 {
91 PrintPutValueError(index, theValue, "PutValue(..) ");
92 }
93 else
94 {
95 dataVector[index] = theValue;
96 }
97 }
98
99 // ---------------------------------------------------------------
100 inline G4PhysicsVectorType G4PhysicsVector::GetType() const
101 {
102 return type;
103 }
104
105 // ---------------------------------------------------------------
106 inline G4bool G4PhysicsVector::GetSpline() const
107 {
108 return useSpline;
109 }
110
111 // ---------------------------------------------------------------
112 inline void G4PhysicsVector::SetVerboseLevel(G4int value)
113 {
114 verboseLevel = value;
115 }
116
117 // ---------------------------------------------------------------
118 inline G4double
119 G4PhysicsVector::FindLinearEnergy(const G4double rand) const
120 {
121 return GetEnergy(rand*dataVector[numberOfNodes - 1]);
122 }
123
124 // ---------------------------------------------------------------
125 inline G4double G4PhysicsVector::Interpolation(const std::size_t idx,
126 const G4double e) const
127 {
128 // perform the interpolation
129 const G4double x1 = binVector[idx];
130 const G4double dl = binVector[idx + 1] - x1;
131
132 const G4double y1 = dataVector[idx];
133 const G4double dy = dataVector[idx + 1] - y1;
134
135 // note: all corner cases of the previous methods are covered and eventually
136 // gives b=0/1 that results in y=y0\y_{N-1} if e<=x[0]/e>=x[N-1] or
137 // y=y_i/y_{i+1} if e<x[i]/e>=x[i+1] due to small numerical errors
138 const G4double b = (e - x1) / dl;
139
140 G4double res = y1 + b * dy;
141
142 if (useSpline) // spline interpolation
143 {
144 const G4double c0 = (2.0 - b) * secDerivative[idx];
145 const G4double c1 = (1.0 + b) * secDerivative[idx + 1];
146 res += (b * (b - 1.0)) * (c0 + c1) * (dl * dl * (1.0/6.0));
147 }
148
149 return res;
150 }
151
152 // ---------------------------------------------------------------
153 inline std::size_t G4PhysicsVector::ComputeLogVectorBin(
154 const G4double loge) const
155 {
156 return static_cast<std::size_t>( std::min( static_cast<G4int>((loge - logemin) * invdBin),
157 static_cast<G4int>(idxmax) ) );
158 }
159
160 // ---------------------------------------------------------------
161 inline std::size_t
162 G4PhysicsVector::LogBin(const G4double e, const G4double loge) const
163 {
164 std::size_t idx =
165 scale[std::min( static_cast<G4int>((loge - lmin1) * iBin1),
166 static_cast<G4int>(imax1) )];
167 for (; idx <= idxmax; ++idx)
168 {
169 if (e >= binVector[idx] && e <= binVector[idx + 1]) { break; }
170 }
171 return idx;
172 }
173
174 // ---------------------------------------------------------------
175 inline std::size_t G4PhysicsVector::BinaryBin(const G4double e) const
176 {
177 // Bin location proposed by K.Genser (FNAL)
178 return std::lower_bound(binVector.cbegin(), binVector.cend(), e) -
179 binVector.cbegin() - 1;
180 }
181
182 // ---------------------------------------------------------------
183 inline std::size_t G4PhysicsVector::GetBin(const G4double e) const
184 {
185 std::size_t bin;
186 switch(type)
187 {
188 case T_G4PhysicsLogVector:
189 bin = ComputeLogVectorBin(G4Log(e));
190 break;
191
192 case T_G4PhysicsLinearVector:
193 bin = static_cast<std::size_t>( std::min( static_cast<G4int>((e - edgeMin) * invdBin),
194 static_cast<G4int>(idxmax) ) );
195 break;
196
197 default:
198 bin = (nLogNodes > 0) ? LogBin(e, G4Log(e)) : BinaryBin(e);
199 }
200 return bin;
201 }
202
203 // ---------------------------------------------------------------
204 inline G4double
205 G4PhysicsVector::Value(const G4double e, std::size_t& idx) const
206 {
207 G4double res;
208 if (idx + 1 < numberOfNodes &&
209 e >= binVector[idx] && e <= binVector[idx+1])
210 {
211 res = Interpolation(idx, e);
212 }
213 else if (e > edgeMin && e < edgeMax)
214 {
215 idx = GetBin(e);
216 res = Interpolation(idx, e);
217 }
218 else if(e <= edgeMin)
219 {
220 res = dataVector[0];
221 idx = 0;
222 }
223 else
224 {
225 res = dataVector[idxmax + 1];
226 idx = idxmax;
227 }
228 return res;
229 }
230
231 // ---------------------------------------------------------------
232 inline G4double G4PhysicsVector::Value(G4double e) const
233 {
234 G4double res;
235 if (e > edgeMin && e < edgeMax)
236 {
237 const std::size_t idx = GetBin(e);
238 res = Interpolation(idx, e);
239 }
240 else if(e <= edgeMin)
241 {
242 res = dataVector[0];
243 }
244 else
245 {
246 res = dataVector[idxmax + 1];
247 }
248 return res;
249 }
250
251 // ---------------------------------------------------------------
252 inline G4double G4PhysicsVector::GetValue(G4double e, G4bool&) const
253 {
254 return Value(e);
255 }
256
257 // ---------------------------------------------------------------
258 inline G4double
259 G4PhysicsVector::LogVectorValue(const G4double e, const G4double loge) const
260 {
261 G4double res;
262 if (e > edgeMin && e < edgeMax)
263 {
264 const std::size_t idx = ComputeLogVectorBin(loge);
265 res = Interpolation(idx, e);
266 }
267 else if (e <= edgeMin)
268 {
269 res = dataVector[0];
270 }
271 else
272 {
273 res = dataVector[idxmax - 1];
274 }
275 return res;
276 }
277
278 // ---------------------------------------------------------------
279 inline G4double
280 G4PhysicsVector::LogFreeVectorValue(const G4double e, const G4double loge) const
281 {
282 G4double res;
283 if (e > edgeMin && e < edgeMax)
284 {
285 const std::size_t idx = LogBin(e, loge);
286 res = Interpolation(idx, e);
287 }
288 else if (e <= edgeMin)
289 {
290 res = dataVector[0];
291 }
292 else
293 {
294 res = dataVector[idxmax + 1];
295 }
296 return res;
297 }
298
299 // ---------------------------------------------------------------
300