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