| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4Physics2DVector class implementation
27 //
28 // Author: Vladimir Ivanchenko, 25.09.2011
29 // -------------------------------------------
30
31 #include <iomanip>
32
33 #include "G4Physics2DVector.hh"
34
35 // -------------------------------------------
36
37 G4Physics2DVector::G4Physics2DVector() { Prepa
38
39 // -------------------------------------------
40
41 G4Physics2DVector::G4Physics2DVector(std::size
42 {
43 if(nx < 2 || ny < 2)
44 {
45 G4ExceptionDescription ed;
46 ed << "G4Physics2DVector is too short: nx=
47 G4Exception("G4Physics2DVector::G4Physics2
48 FatalException, ed, "Both leng
49 }
50 numberOfXNodes = nx;
51 numberOfYNodes = ny;
52 PrepareVectors();
53 }
54
55 // -------------------------------------------
56
57 G4Physics2DVector::~G4Physics2DVector() { Clea
58
59 // -------------------------------------------
60
61 G4Physics2DVector::G4Physics2DVector(const G4P
62 {
63 type = right.type;
64
65 numberOfXNodes = right.numberOfXNodes;
66 numberOfYNodes = right.numberOfYNodes;
67
68 verboseLevel = right.verboseLevel;
69 useBicubic = right.useBicubic;
70
71 xVector = right.xVector;
72 yVector = right.yVector;
73
74 PrepareVectors();
75 CopyData(right);
76 }
77
78 // -------------------------------------------
79
80 G4Physics2DVector& G4Physics2DVector::operator
81 {
82 if(&right == this)
83 {
84 return *this;
85 }
86 ClearVectors();
87
88 type = right.type;
89
90 numberOfXNodes = right.numberOfXNodes;
91 numberOfYNodes = right.numberOfYNodes;
92
93 verboseLevel = right.verboseLevel;
94 useBicubic = right.useBicubic;
95
96 PrepareVectors();
97 CopyData(right);
98
99 return *this;
100 }
101
102 // -------------------------------------------
103
104 void G4Physics2DVector::PrepareVectors()
105 {
106 xVector.resize(numberOfXNodes, 0.);
107 yVector.resize(numberOfYNodes, 0.);
108 value.resize(numberOfYNodes);
109 for(std::size_t j = 0; j < numberOfYNodes; +
110 {
111 value[j] = new G4PV2DDataVector(numberOfXN
112 }
113 }
114
115 // -------------------------------------------
116
117 void G4Physics2DVector::ClearVectors()
118 {
119 for(std::size_t j = 0; j < numberOfYNodes; +
120 {
121 delete value[j];
122 }
123 }
124
125 // -------------------------------------------
126
127 void G4Physics2DVector::CopyData(const G4Physi
128 {
129 for(std::size_t i = 0; i < numberOfXNodes; +
130 {
131 xVector[i] = right.xVector[i];
132 }
133 for(std::size_t j = 0; j < numberOfYNodes; +
134 {
135 yVector[j] = right.yVector[j];
136 G4PV2DDataVector* v0 = right.value[j];
137 for(std::size_t i = 0; i < numberOfXNodes;
138 {
139 PutValue(i, j, (*v0)[i]);
140 }
141 }
142 }
143
144 // -------------------------------------------
145
146 G4double G4Physics2DVector::Value(G4double xx,
147 std::size_t&
148 {
149 // no interpolation outside the table
150 const G4double x =
151 std::min(std::max(xx, xVector[0]), xVector
152 const G4double y =
153 std::min(std::max(yy, yVector[0]), yVector
154
155 // find bins
156 idx = FindBinLocationX(x, idx);
157 idy = FindBinLocationY(y, idy);
158
159 // interpolate
160 if(useBicubic)
161 {
162 return BicubicInterpolation(x, y, idx, idy
163 }
164
165 const G4double x1 = xVector[idx];
166 const G4double x2 = xVector[idx + 1];
167 const G4double y1 = yVector[idy];
168 const G4double y2 = yVector[idy + 1];
169 const G4double v11 = GetValue(idx, idy);
170 const G4double v12 = GetValue(idx + 1, idy);
171 const G4double v21 = GetValue(idx, idy + 1);
172 const G4double v22 = GetValue(idx + 1, idy +
173 return ((y2 - y) * (v11 * (x2 - x) + v12 * (
174 ((y - y1) * (v21 * (x2 - x) + v22 *
175 ((x2 - x1) * (y2 - y1));
176
177 }
178
179 // -------------------------------------------
180
181 G4double G4Physics2DVector::BicubicInterpolati
182
183
184
185 {
186 // Bicubic interpolation according to
187 // 1. H.M. Antia, "Numerical Methods for Sci
188 // MGH, 1991.
189 // 2. W.H. Press et al., "Numerical recipes.
190 // Computing", Cambridge University Press
191 const G4double x1 = xVector[idx];
192 const G4double x2 = xVector[idx + 1];
193 const G4double y1 = yVector[idy];
194 const G4double y2 = yVector[idy + 1];
195 const G4double f1 = GetValue(idx, idy);
196 const G4double f2 = GetValue(idx + 1, idy);
197 const G4double f3 = GetValue(idx + 1, idy +
198 const G4double f4 = GetValue(idx, idy + 1);
199
200 const G4double dx = x2 - x1;
201 const G4double dy = y2 - y1;
202
203 const G4double h1 = (x - x1) / dx;
204 const G4double h2 = (y - y1) / dy;
205
206 const G4double h12 = h1 * h1;
207 const G4double h13 = h12 * h1;
208 const G4double h22 = h2 * h2;
209 const G4double h23 = h22 * h2;
210
211 // Three derivatives at each of four points
212 // subregion are computed by numerical cente
213 // the functional values already tabulated o
214
215 const G4double f1x = DerivativeX(idx, idy, d
216 const G4double f2x = DerivativeX(idx + 1, id
217 const G4double f3x = DerivativeX(idx + 1, id
218 const G4double f4x = DerivativeX(idx, idy +
219
220 const G4double f1y = DerivativeY(idx, idy, d
221 const G4double f2y = DerivativeY(idx + 1, id
222 const G4double f3y = DerivativeY(idx + 1, id
223 const G4double f4y = DerivativeY(idx, idy +
224
225 const G4double dxy = dx * dy;
226 const G4double f1xy = DerivativeXY(idx, idy,
227 const G4double f2xy = DerivativeXY(idx + 1,
228 const G4double f3xy = DerivativeXY(idx + 1,
229 const G4double f4xy = DerivativeXY(idx, idy
230
231 return f1 + f1y * h2 + (3 * (f4 - f1) - 2 *
232 (2 * (f1 - f4) + f1y + f4y) * h23 + f
233 (3 * (f4x - f1x) - 2 * f1xy - f4xy) *
234 (2 * (f1x - f4x) + f1xy + f4xy) * h1
235 (3 * (f2 - f1) - 2 * f1x - f2x) * h12
236 (3 * f2y - 3 * f1y - 2 * f1xy - f2xy)
237 (9 * (f1 - f2 + f3 - f4) + 6 * f1x +
238 6 * f1y - 6 * f2y - 3 * f3y + 3 * f4
239 2 * f4xy) *
240 h12 * h22 +
241 (6 * (-f1 + f2 - f3 + f4) - 4 * f1x -
242 3 * f1y + 3 * f2y + 3 * f3y - 3 * f4
243 2 * f4xy) *
244 h12 * h23 +
245 (2 * (f1 - f2) + f1x + f2x) * h13 +
246 (2 * (f1y - f2y) + f1xy + f2xy) * h13
247 (6 * (-f1 + f2 - f3 + f4) + 3 * (-f1x
248 4 * f2y + 2 * f3y - 2 * f4y - 2 * f1
249 h13 * h22 +
250 (4 * (f1 - f2 + f3 - f4) + 2 * (f1x +
251 2 * (f1y - f2y - f3y + f4y) + f1xy +
252 h13 * h23;
253 }
254
255 // -------------------------------------------
256
257 void G4Physics2DVector::PutVectors(const std::
258 const std::
259 {
260 ClearVectors();
261 std::size_t nx = vecX.size();
262 std::size_t ny = vecY.size();
263 if(nx < 2 || ny < 2)
264 {
265 G4ExceptionDescription ed;
266 ed << "G4Physics2DVector is too short: nx=
267 G4Exception("G4Physics2DVector::PutVectors
268 "Both lengths should be above
269 }
270
271 numberOfXNodes = nx;
272 numberOfYNodes = ny;
273 PrepareVectors();
274 for(std::size_t i = 0; i < nx; ++i)
275 {
276 xVector[i] = vecX[i];
277 }
278 for(std::size_t j = 0; j < ny; ++j)
279 {
280 yVector[j] = vecY[j];
281 }
282 }
283
284 // -------------------------------------------
285
286 void G4Physics2DVector::Store(std::ofstream& o
287 {
288 // binning
289 G4long prec = out.precision();
290 out << G4int(type) << " " << numberOfXNodes
291 << G4endl;
292 out << std::setprecision(8);
293
294 // contents
295 for(std::size_t i = 0; i < numberOfXNodes -
296 {
297 out << xVector[i] << " ";
298 }
299 out << xVector[numberOfXNodes - 1] << G4endl
300 for(std::size_t j = 0; j < numberOfYNodes -
301 {
302 out << yVector[j] << " ";
303 }
304 out << yVector[numberOfYNodes - 1] << G4endl
305 for(std::size_t j = 0; j < numberOfYNodes; +
306 {
307 for(std::size_t i = 0; i < numberOfXNodes
308 {
309 out << GetValue(i, j) << " ";
310 }
311 out << GetValue(numberOfXNodes - 1, j) <<
312 }
313 out.precision(prec);
314 out.close();
315 }
316
317 // -------------------------------------------
318
319 G4bool G4Physics2DVector::Retrieve(std::ifstre
320 {
321 // initialisation
322 ClearVectors();
323
324 // binning
325 G4int k, nx, ny;
326 in >> k >> nx >> ny;
327 if(in.fail() || 2 > nx || 2 > ny || nx >= IN
328 {
329 return false;
330 }
331 numberOfXNodes = nx;
332 numberOfYNodes = ny;
333 PrepareVectors();
334 type = G4PhysicsVectorType(k);
335
336 // contents
337 G4double val;
338 for(G4int i = 0; i < nx; ++i)
339 {
340 in >> xVector[i];
341 if(in.fail())
342 {
343 return false;
344 }
345 }
346 for(G4int j = 0; j < ny; ++j)
347 {
348 in >> yVector[j];
349 if(in.fail())
350 {
351 return false;
352 }
353 }
354 for(G4int j = 0; j < ny; ++j)
355 {
356 for(G4int i = 0; i < nx; ++i)
357 {
358 in >> val;
359 if(in.fail())
360 {
361 return false;
362 }
363 PutValue(i, j, val);
364 }
365 }
366 in.close();
367 return true;
368 }
369
370 // -------------------------------------------
371
372 void G4Physics2DVector::ScaleVector(G4double f
373 {
374 G4double val;
375 for(std::size_t j = 0; j < numberOfYNodes; +
376 {
377 for(std::size_t i = 0; i < numberOfXNodes;
378 {
379 val = GetValue(i, j) * factor;
380 PutValue(i, j, val);
381 }
382 }
383 }
384
385 // -------------------------------------------
386
387 G4double G4Physics2DVector::FindLinearX(G4doub
388 std::s
389 {
390 G4double y = std::min(std::max(yy, yVector[0
391
392 // find bins
393 idy = FindBinLocationY(y, idy);
394
395 G4double x1 = InterpolateLinearX(*(value[id
396 G4double x2 = InterpolateLinearX(*(value[id
397 G4double res = x1;
398 G4double del = yVector[idy + 1] - yVector[id
399 if(del != 0.0)
400 {
401 res += (x2 - x1) * (y - yVector[idy]) / de
402 }
403 return res;
404 }
405
406 // -------------------------------------------
407
408 G4double G4Physics2DVector::InterpolateLinearX
409
410 {
411 std::size_t nn = v.size();
412 if(1 >= nn)
413 {
414 return 0.0;
415 }
416 std::size_t n1 = 0;
417 std::size_t n2 = nn / 2;
418 std::size_t n3 = nn - 1;
419 G4double y = rand * v[n3];
420 while(n1 + 1 != n3)
421 {
422 if(y > v[n2])
423 {
424 n1 = n2;
425 }
426 else
427 {
428 n3 = n2;
429 }
430 n2 = (n3 + n1 + 1) / 2;
431 }
432 G4double res = xVector[n1];
433 G4double del = v[n3] - v[n1];
434 if(del > 0.0)
435 {
436 res += (y - v[n1]) * (xVector[n3] - res) /
437 }
438 return res;
439 }
440
441 // -------------------------------------------
442