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 165 const G4double x1 = xVector[idx]; 165 const G4double x1 = xVector[idx]; 166 const G4double x2 = xVector[idx + 1]; 166 const G4double x2 = xVector[idx + 1]; 167 const G4double y1 = yVector[idy]; 167 const G4double y1 = yVector[idy]; 168 const G4double y2 = yVector[idy + 1]; 168 const G4double y2 = yVector[idy + 1]; 169 const G4double v11 = GetValue(idx, idy); 169 const G4double v11 = GetValue(idx, idy); 170 const G4double v12 = GetValue(idx + 1, idy); 170 const G4double v12 = GetValue(idx + 1, idy); 171 const G4double v21 = GetValue(idx, idy + 1); 171 const G4double v21 = GetValue(idx, idy + 1); 172 const G4double v22 = GetValue(idx + 1, idy + 172 const G4double v22 = GetValue(idx + 1, idy + 1); 173 return ((y2 - y) * (v11 * (x2 - x) + v12 * ( 173 return ((y2 - y) * (v11 * (x2 - x) + v12 * (x - x1)) + 174 ((y - y1) * (v21 * (x2 - x) + v22 * 174 ((y - y1) * (v21 * (x2 - x) + v22 * (x - x1)))) / 175 ((x2 - x1) * (y2 - y1)); 175 ((x2 - x1) * (y2 - y1)); 176 176 177 } 177 } 178 178 179 // ------------------------------------------- 179 // -------------------------------------------------------------- 180 180 181 G4double G4Physics2DVector::BicubicInterpolati 181 G4double G4Physics2DVector::BicubicInterpolation(const G4double x, 182 182 const G4double y, 183 183 const std::size_t idx, 184 184 const std::size_t idy) const 185 { 185 { 186 // Bicubic interpolation according to 186 // Bicubic interpolation according to 187 // 1. H.M. Antia, "Numerical Methods for Sci 187 // 1. H.M. Antia, "Numerical Methods for Scientists and Engineers", 188 // MGH, 1991. 188 // MGH, 1991. 189 // 2. W.H. Press et al., "Numerical recipes. 189 // 2. W.H. Press et al., "Numerical recipes. The Art of Scientific 190 // Computing", Cambridge University Press 190 // Computing", Cambridge University Press, 2007. 191 const G4double x1 = xVector[idx]; 191 const G4double x1 = xVector[idx]; 192 const G4double x2 = xVector[idx + 1]; 192 const G4double x2 = xVector[idx + 1]; 193 const G4double y1 = yVector[idy]; 193 const G4double y1 = yVector[idy]; 194 const G4double y2 = yVector[idy + 1]; 194 const G4double y2 = yVector[idy + 1]; 195 const G4double f1 = GetValue(idx, idy); 195 const G4double f1 = GetValue(idx, idy); 196 const G4double f2 = GetValue(idx + 1, idy); 196 const G4double f2 = GetValue(idx + 1, idy); 197 const G4double f3 = GetValue(idx + 1, idy + 197 const G4double f3 = GetValue(idx + 1, idy + 1); 198 const G4double f4 = GetValue(idx, idy + 1); 198 const G4double f4 = GetValue(idx, idy + 1); 199 199 200 const G4double dx = x2 - x1; 200 const G4double dx = x2 - x1; 201 const G4double dy = y2 - y1; 201 const G4double dy = y2 - y1; 202 202 203 const G4double h1 = (x - x1) / dx; 203 const G4double h1 = (x - x1) / dx; 204 const G4double h2 = (y - y1) / dy; 204 const G4double h2 = (y - y1) / dy; 205 205 206 const G4double h12 = h1 * h1; 206 const G4double h12 = h1 * h1; 207 const G4double h13 = h12 * h1; 207 const G4double h13 = h12 * h1; 208 const G4double h22 = h2 * h2; 208 const G4double h22 = h2 * h2; 209 const G4double h23 = h22 * h2; 209 const G4double h23 = h22 * h2; 210 210 211 // Three derivatives at each of four points 211 // Three derivatives at each of four points (1-4) defining the 212 // subregion are computed by numerical cente 212 // subregion are computed by numerical centered differencing from 213 // the functional values already tabulated o 213 // the functional values already tabulated on the grid. 214 214 215 const G4double f1x = DerivativeX(idx, idy, d 215 const G4double f1x = DerivativeX(idx, idy, dx); 216 const G4double f2x = DerivativeX(idx + 1, id 216 const G4double f2x = DerivativeX(idx + 1, idy, dx); 217 const G4double f3x = DerivativeX(idx + 1, id 217 const G4double f3x = DerivativeX(idx + 1, idy + 1, dx); 218 const G4double f4x = DerivativeX(idx, idy + 218 const G4double f4x = DerivativeX(idx, idy + 1, dx); 219 219 220 const G4double f1y = DerivativeY(idx, idy, d 220 const G4double f1y = DerivativeY(idx, idy, dy); 221 const G4double f2y = DerivativeY(idx + 1, id 221 const G4double f2y = DerivativeY(idx + 1, idy, dy); 222 const G4double f3y = DerivativeY(idx + 1, id 222 const G4double f3y = DerivativeY(idx + 1, idy + 1, dy); 223 const G4double f4y = DerivativeY(idx, idy + 223 const G4double f4y = DerivativeY(idx, idy + 1, dy); 224 224 225 const G4double dxy = dx * dy; 225 const G4double dxy = dx * dy; 226 const G4double f1xy = DerivativeXY(idx, idy, 226 const G4double f1xy = DerivativeXY(idx, idy, dxy); 227 const G4double f2xy = DerivativeXY(idx + 1, 227 const G4double f2xy = DerivativeXY(idx + 1, idy, dxy); 228 const G4double f3xy = DerivativeXY(idx + 1, 228 const G4double f3xy = DerivativeXY(idx + 1, idy + 1, dxy); 229 const G4double f4xy = DerivativeXY(idx, idy 229 const G4double f4xy = DerivativeXY(idx, idy + 1, dxy); 230 230 231 return f1 + f1y * h2 + (3 * (f4 - f1) - 2 * 231 return f1 + f1y * h2 + (3 * (f4 - f1) - 2 * f1y - f4y) * h22 + 232 (2 * (f1 - f4) + f1y + f4y) * h23 + f 232 (2 * (f1 - f4) + f1y + f4y) * h23 + f1x * h1 + f1xy * h1 * h2 + 233 (3 * (f4x - f1x) - 2 * f1xy - f4xy) * 233 (3 * (f4x - f1x) - 2 * f1xy - f4xy) * h1 * h22 + 234 (2 * (f1x - f4x) + f1xy + f4xy) * h1 234 (2 * (f1x - f4x) + f1xy + f4xy) * h1 * h23 + 235 (3 * (f2 - f1) - 2 * f1x - f2x) * h12 235 (3 * (f2 - f1) - 2 * f1x - f2x) * h12 + 236 (3 * f2y - 3 * f1y - 2 * f1xy - f2xy) 236 (3 * f2y - 3 * f1y - 2 * f1xy - f2xy) * h12 * h2 + 237 (9 * (f1 - f2 + f3 - f4) + 6 * f1x + 237 (9 * (f1 - f2 + f3 - f4) + 6 * f1x + 3 * f2x - 3 * f3x - 6 * f4x + 238 6 * f1y - 6 * f2y - 3 * f3y + 3 * f4 238 6 * f1y - 6 * f2y - 3 * f3y + 3 * f4y + 4 * f1xy + 2 * f2xy + f3xy + 239 2 * f4xy) * 239 2 * f4xy) * 240 h12 * h22 + 240 h12 * h22 + 241 (6 * (-f1 + f2 - f3 + f4) - 4 * f1x - 241 (6 * (-f1 + f2 - f3 + f4) - 4 * f1x - 2 * f2x + 2 * f3x + 4 * f4x - 242 3 * f1y + 3 * f2y + 3 * f3y - 3 * f4 242 3 * f1y + 3 * f2y + 3 * f3y - 3 * f4y - 2 * f1xy - f2xy - f3xy - 243 2 * f4xy) * 243 2 * f4xy) * 244 h12 * h23 + 244 h12 * h23 + 245 (2 * (f1 - f2) + f1x + f2x) * h13 + 245 (2 * (f1 - f2) + f1x + f2x) * h13 + 246 (2 * (f1y - f2y) + f1xy + f2xy) * h13 246 (2 * (f1y - f2y) + f1xy + f2xy) * h13 * h2 + 247 (6 * (-f1 + f2 - f3 + f4) + 3 * (-f1x 247 (6 * (-f1 + f2 - f3 + f4) + 3 * (-f1x - f2x + f3x + f4x) - 4 * f1y + 248 4 * f2y + 2 * f3y - 2 * f4y - 2 * f1 248 4 * f2y + 2 * f3y - 2 * f4y - 2 * f1xy - 2 * f2xy - f3xy - f4xy) * 249 h13 * h22 + 249 h13 * h22 + 250 (4 * (f1 - f2 + f3 - f4) + 2 * (f1x + 250 (4 * (f1 - f2 + f3 - f4) + 2 * (f1x + f2x - f3x - f4x) + 251 2 * (f1y - f2y - f3y + f4y) + f1xy + 251 2 * (f1y - f2y - f3y + f4y) + f1xy + f2xy + f3xy + f4xy) * 252 h13 * h23; 252 h13 * h23; 253 } 253 } 254 254 255 // ------------------------------------------- 255 // -------------------------------------------------------------- 256 256 257 void G4Physics2DVector::PutVectors(const std:: 257 void G4Physics2DVector::PutVectors(const std::vector<G4double>& vecX, 258 const std:: 258 const std::vector<G4double>& vecY) 259 { 259 { 260 ClearVectors(); 260 ClearVectors(); 261 std::size_t nx = vecX.size(); 261 std::size_t nx = vecX.size(); 262 std::size_t ny = vecY.size(); 262 std::size_t ny = vecY.size(); 263 if(nx < 2 || ny < 2) 263 if(nx < 2 || ny < 2) 264 { 264 { 265 G4ExceptionDescription ed; 265 G4ExceptionDescription ed; 266 ed << "G4Physics2DVector is too short: nx= 266 ed << "G4Physics2DVector is too short: nx= " << nx << " ny= " << ny; 267 G4Exception("G4Physics2DVector::PutVectors 267 G4Exception("G4Physics2DVector::PutVectors()", "glob03", FatalException, ed, 268 "Both lengths should be above 268 "Both lengths should be above 1"); 269 } 269 } 270 270 271 numberOfXNodes = nx; 271 numberOfXNodes = nx; 272 numberOfYNodes = ny; 272 numberOfYNodes = ny; 273 PrepareVectors(); 273 PrepareVectors(); 274 for(std::size_t i = 0; i < nx; ++i) 274 for(std::size_t i = 0; i < nx; ++i) 275 { 275 { 276 xVector[i] = vecX[i]; 276 xVector[i] = vecX[i]; 277 } 277 } 278 for(std::size_t j = 0; j < ny; ++j) 278 for(std::size_t j = 0; j < ny; ++j) 279 { 279 { 280 yVector[j] = vecY[j]; 280 yVector[j] = vecY[j]; 281 } 281 } 282 } 282 } 283 283 284 // ------------------------------------------- 284 // -------------------------------------------------------------- 285 285 286 void G4Physics2DVector::Store(std::ofstream& o 286 void G4Physics2DVector::Store(std::ofstream& out) const 287 { 287 { 288 // binning 288 // binning 289 G4long prec = out.precision(); 289 G4long prec = out.precision(); 290 out << G4int(type) << " " << numberOfXNodes 290 out << G4int(type) << " " << numberOfXNodes << " " << numberOfYNodes 291 << G4endl; 291 << G4endl; 292 out << std::setprecision(8); 292 out << std::setprecision(8); 293 293 294 // contents 294 // contents 295 for(std::size_t i = 0; i < numberOfXNodes - 295 for(std::size_t i = 0; i < numberOfXNodes - 1; ++i) 296 { 296 { 297 out << xVector[i] << " "; 297 out << xVector[i] << " "; 298 } 298 } 299 out << xVector[numberOfXNodes - 1] << G4endl 299 out << xVector[numberOfXNodes - 1] << G4endl; 300 for(std::size_t j = 0; j < numberOfYNodes - 300 for(std::size_t j = 0; j < numberOfYNodes - 1; ++j) 301 { 301 { 302 out << yVector[j] << " "; 302 out << yVector[j] << " "; 303 } 303 } 304 out << yVector[numberOfYNodes - 1] << G4endl 304 out << yVector[numberOfYNodes - 1] << G4endl; 305 for(std::size_t j = 0; j < numberOfYNodes; + 305 for(std::size_t j = 0; j < numberOfYNodes; ++j) 306 { 306 { 307 for(std::size_t i = 0; i < numberOfXNodes 307 for(std::size_t i = 0; i < numberOfXNodes - 1; ++i) 308 { 308 { 309 out << GetValue(i, j) << " "; 309 out << GetValue(i, j) << " "; 310 } 310 } 311 out << GetValue(numberOfXNodes - 1, j) << 311 out << GetValue(numberOfXNodes - 1, j) << G4endl; 312 } 312 } 313 out.precision(prec); 313 out.precision(prec); 314 out.close(); 314 out.close(); 315 } 315 } 316 316 317 // ------------------------------------------- 317 // -------------------------------------------------------------- 318 318 319 G4bool G4Physics2DVector::Retrieve(std::ifstre 319 G4bool G4Physics2DVector::Retrieve(std::ifstream& in) 320 { 320 { 321 // initialisation 321 // initialisation 322 ClearVectors(); 322 ClearVectors(); 323 323 324 // binning 324 // binning 325 G4int k, nx, ny; 325 G4int k, nx, ny; 326 in >> k >> nx >> ny; 326 in >> k >> nx >> ny; 327 if(in.fail() || 2 > nx || 2 > ny || nx >= IN 327 if(in.fail() || 2 > nx || 2 > ny || nx >= INT_MAX || ny >= INT_MAX) 328 { 328 { 329 return false; 329 return false; 330 } 330 } 331 numberOfXNodes = nx; 331 numberOfXNodes = nx; 332 numberOfYNodes = ny; 332 numberOfYNodes = ny; 333 PrepareVectors(); 333 PrepareVectors(); 334 type = G4PhysicsVectorType(k); 334 type = G4PhysicsVectorType(k); 335 335 336 // contents 336 // contents 337 G4double val; 337 G4double val; 338 for(G4int i = 0; i < nx; ++i) 338 for(G4int i = 0; i < nx; ++i) 339 { 339 { 340 in >> xVector[i]; 340 in >> xVector[i]; 341 if(in.fail()) 341 if(in.fail()) 342 { 342 { 343 return false; 343 return false; 344 } 344 } 345 } 345 } 346 for(G4int j = 0; j < ny; ++j) 346 for(G4int j = 0; j < ny; ++j) 347 { 347 { 348 in >> yVector[j]; 348 in >> yVector[j]; 349 if(in.fail()) 349 if(in.fail()) 350 { 350 { 351 return false; 351 return false; 352 } 352 } 353 } 353 } 354 for(G4int j = 0; j < ny; ++j) 354 for(G4int j = 0; j < ny; ++j) 355 { 355 { 356 for(G4int i = 0; i < nx; ++i) 356 for(G4int i = 0; i < nx; ++i) 357 { 357 { 358 in >> val; 358 in >> val; 359 if(in.fail()) 359 if(in.fail()) 360 { 360 { 361 return false; 361 return false; 362 } 362 } 363 PutValue(i, j, val); 363 PutValue(i, j, val); 364 } 364 } 365 } 365 } 366 in.close(); 366 in.close(); 367 return true; 367 return true; 368 } 368 } 369 369 370 // ------------------------------------------- 370 // -------------------------------------------------------------- 371 371 372 void G4Physics2DVector::ScaleVector(G4double f 372 void G4Physics2DVector::ScaleVector(G4double factor) 373 { 373 { 374 G4double val; 374 G4double val; 375 for(std::size_t j = 0; j < numberOfYNodes; + 375 for(std::size_t j = 0; j < numberOfYNodes; ++j) 376 { 376 { 377 for(std::size_t i = 0; i < numberOfXNodes; 377 for(std::size_t i = 0; i < numberOfXNodes; ++i) 378 { 378 { 379 val = GetValue(i, j) * factor; 379 val = GetValue(i, j) * factor; 380 PutValue(i, j, val); 380 PutValue(i, j, val); 381 } 381 } 382 } 382 } 383 } 383 } 384 384 385 // ------------------------------------------- 385 // -------------------------------------------------------------- 386 386 387 G4double G4Physics2DVector::FindLinearX(G4doub 387 G4double G4Physics2DVector::FindLinearX(G4double rand, G4double yy, 388 std::s 388 std::size_t& idy) const 389 { 389 { 390 G4double y = std::min(std::max(yy, yVector[0 390 G4double y = std::min(std::max(yy, yVector[0]), yVector[numberOfYNodes - 1]); 391 391 392 // find bins 392 // find bins 393 idy = FindBinLocationY(y, idy); 393 idy = FindBinLocationY(y, idy); 394 394 395 G4double x1 = InterpolateLinearX(*(value[id 395 G4double x1 = InterpolateLinearX(*(value[idy]), rand); 396 G4double x2 = InterpolateLinearX(*(value[id 396 G4double x2 = InterpolateLinearX(*(value[idy + 1]), rand); 397 G4double res = x1; 397 G4double res = x1; 398 G4double del = yVector[idy + 1] - yVector[id 398 G4double del = yVector[idy + 1] - yVector[idy]; 399 if(del != 0.0) 399 if(del != 0.0) 400 { 400 { 401 res += (x2 - x1) * (y - yVector[idy]) / de 401 res += (x2 - x1) * (y - yVector[idy]) / del; 402 } 402 } 403 return res; 403 return res; 404 } 404 } 405 405 406 // ------------------------------------------- 406 // -------------------------------------------------------------- 407 407 408 G4double G4Physics2DVector::InterpolateLinearX 408 G4double G4Physics2DVector::InterpolateLinearX(G4PV2DDataVector& v, 409 409 G4double rand) const 410 { 410 { 411 std::size_t nn = v.size(); 411 std::size_t nn = v.size(); 412 if(1 >= nn) 412 if(1 >= nn) 413 { 413 { 414 return 0.0; 414 return 0.0; 415 } 415 } 416 std::size_t n1 = 0; 416 std::size_t n1 = 0; 417 std::size_t n2 = nn / 2; 417 std::size_t n2 = nn / 2; 418 std::size_t n3 = nn - 1; 418 std::size_t n3 = nn - 1; 419 G4double y = rand * v[n3]; 419 G4double y = rand * v[n3]; 420 while(n1 + 1 != n3) 420 while(n1 + 1 != n3) 421 { 421 { 422 if(y > v[n2]) 422 if(y > v[n2]) 423 { 423 { 424 n1 = n2; 424 n1 = n2; 425 } 425 } 426 else 426 else 427 { 427 { 428 n3 = n2; 428 n3 = n2; 429 } 429 } 430 n2 = (n3 + n1 + 1) / 2; 430 n2 = (n3 + n1 + 1) / 2; 431 } 431 } 432 G4double res = xVector[n1]; 432 G4double res = xVector[n1]; 433 G4double del = v[n3] - v[n1]; 433 G4double del = v[n3] - v[n1]; 434 if(del > 0.0) 434 if(del > 0.0) 435 { 435 { 436 res += (y - v[n1]) * (xVector[n3] - res) / 436 res += (y - v[n1]) * (xVector[n3] - res) / del; 437 } 437 } 438 return res; 438 return res; 439 } 439 } 440 440 441 // ------------------------------------------- 441 // -------------------------------------------------------------- 442 442