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