00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "G4Physics2DVector.hh"
00043 #include <iomanip>
00044
00045
00046
00047 G4Physics2DVector::G4Physics2DVector()
00048 : type(T_G4PhysicsFreeVector),
00049 numberOfXNodes(0), numberOfYNodes(0),
00050 verboseLevel(0), useBicubic(false)
00051 {
00052 cache = new G4Physics2DVectorCache();
00053 }
00054
00055
00056
00057 G4Physics2DVector::G4Physics2DVector(size_t nx, size_t ny)
00058 : type(T_G4PhysicsFreeVector),
00059 numberOfXNodes(nx), numberOfYNodes(ny),
00060 verboseLevel(0), useBicubic(false)
00061 {
00062 cache = new G4Physics2DVectorCache();
00063 PrepareVectors();
00064 }
00065
00066
00067
00068 G4Physics2DVector::~G4Physics2DVector()
00069 {
00070 delete cache;
00071 ClearVectors();
00072 }
00073
00074
00075
00076 G4Physics2DVector::G4Physics2DVector(const G4Physics2DVector& right)
00077 {
00078 type = right.type;
00079
00080 numberOfXNodes = right.numberOfXNodes;
00081 numberOfYNodes = right.numberOfYNodes;
00082
00083 verboseLevel = right.verboseLevel;
00084 useBicubic = right.useBicubic;
00085
00086 xVector = right.xVector;
00087 yVector = right.yVector;
00088
00089 cache = new G4Physics2DVectorCache();
00090 PrepareVectors();
00091 CopyData(right);
00092 }
00093
00094
00095
00096 G4Physics2DVector& G4Physics2DVector::operator=(const G4Physics2DVector& right)
00097 {
00098 if (&right==this) { return *this; }
00099 ClearVectors();
00100
00101 type = right.type;
00102
00103 numberOfXNodes = right.numberOfXNodes;
00104 numberOfYNodes = right.numberOfYNodes;
00105
00106 verboseLevel = right.verboseLevel;
00107 useBicubic = right.useBicubic;
00108
00109 cache->Clear();
00110 PrepareVectors();
00111 CopyData(right);
00112
00113 return *this;
00114 }
00115
00116
00117
00118 void G4Physics2DVector::PrepareVectors()
00119 {
00120 xVector.resize(numberOfXNodes,0.);
00121 yVector.resize(numberOfYNodes,0.);
00122 value.resize(numberOfYNodes,0);
00123 for(size_t j=0; j<numberOfYNodes; ++j) {
00124 G4PV2DDataVector* v = new G4PV2DDataVector();
00125 v->resize(numberOfXNodes,0.);
00126 value[j] = v;
00127 }
00128 }
00129
00130
00131
00132 void G4Physics2DVector::ClearVectors()
00133 {
00134 for(size_t j=0; j<numberOfYNodes; ++j) {
00135 delete value[j];
00136 }
00137 }
00138
00139
00140
00141 void G4Physics2DVector::CopyData(const G4Physics2DVector &right)
00142 {
00143 for(size_t i=0; i<numberOfXNodes; ++i) {
00144 xVector[i] = right.xVector[i];
00145 }
00146 for(size_t j=0; j<numberOfYNodes; ++j) {
00147 yVector[j] = right.yVector[j];
00148 G4PV2DDataVector* v0 = right.value[j];
00149 for(size_t i=0; i<numberOfXNodes; ++i) {
00150 PutValue(i,j,(*v0)[i]);
00151 }
00152 }
00153 }
00154
00155
00156
00157 void G4Physics2DVector::ComputeValue(G4double xx, G4double yy)
00158 {
00159 if(xx != cache->lastBinX) {
00160 if(xx <= xVector[0]) {
00161 cache->lastX = xVector[0];
00162 cache->lastBinX = 0;
00163 } else if(xx >= xVector[numberOfXNodes-1]) {
00164 cache->lastX = xVector[numberOfXNodes-1];
00165 cache->lastBinX = numberOfXNodes-2;
00166 } else {
00167 cache->lastX = xx;
00168 FindBinLocationX(xx);
00169 }
00170 }
00171 if(yy != cache->lastBinY) {
00172 if(yy <= yVector[0]) {
00173 cache->lastY = yVector[0];
00174 cache->lastBinY = 0;
00175 } else if(yy >= yVector[numberOfYNodes-1]) {
00176 cache->lastY = yVector[numberOfYNodes-1];
00177 cache->lastBinY = numberOfYNodes-2;
00178 } else {
00179 cache->lastY = yy;
00180 FindBinLocationY(yy);
00181 }
00182 }
00183 size_t idx = cache->lastBinX;
00184 size_t idy = cache->lastBinY;
00185 if(useBicubic) {
00186 BicubicInterpolation(idx, idy);
00187 } else {
00188 G4double x1 = xVector[idx];
00189 G4double x2 = xVector[idx+1];
00190 G4double y1 = yVector[idy];
00191 G4double y2 = yVector[idy+1];
00192 G4double x = cache->lastX;
00193 G4double y = cache->lastY;
00194 G4double v11= GetValue(idx, idy);
00195 G4double v12= GetValue(idx+1, idy);
00196 G4double v21= GetValue(idx, idy+1);
00197 G4double v22= GetValue(idx+1, idy+1);
00198 cache->lastValue =
00199 ((y2 - y)*(v11*(x2 - x) + v12*(x - x1)) +
00200 ((y - y1)*(v21*(x2 - x) + v22*(x - x1))))/((x2 - x1)*(y2 - y1));
00201 }
00202 }
00203
00204
00205
00206 void G4Physics2DVector::BicubicInterpolation(size_t idx, size_t idy)
00207 {
00208
00209
00210
00211
00212
00213 G4double x1 = xVector[idx];
00214 G4double x2 = xVector[idx+1];
00215 G4double y1 = yVector[idy];
00216 G4double y2 = yVector[idy+1];
00217 G4double x = cache->lastX;
00218 G4double y = cache->lastY;
00219 G4double f1 = GetValue(idx, idy);
00220 G4double f2 = GetValue(idx+1, idy);
00221 G4double f3 = GetValue(idx+1, idy+1);
00222 G4double f4 = GetValue(idx, idy+1);
00223
00224 G4double dx = x2 - x1;
00225 G4double dy = y2 - y1;
00226
00227 G4double h1 = (x - x1)/dx;
00228 G4double h2 = (y - y1)/dy;
00229
00230 G4double h12 = h1*h1;
00231 G4double h13 = h12*h1;
00232 G4double h22 = h2*h2;
00233 G4double h23 = h22*h2;
00234
00235
00236
00237
00238
00239 G4double f1x = DerivativeX(idx, idy, dx);
00240 G4double f2x = DerivativeX(idx+1, idy, dx);
00241 G4double f3x = DerivativeX(idx+1, idy+1, dx);
00242 G4double f4x = DerivativeX(idx, idy+1, dx);
00243
00244 G4double f1y = DerivativeY(idx, idy, dy);
00245 G4double f2y = DerivativeY(idx+1, idy, dy);
00246 G4double f3y = DerivativeY(idx+1, idy+1, dy);
00247 G4double f4y = DerivativeY(idx, idy+1, dy);
00248
00249 G4double dxy = dx*dy;
00250 G4double f1xy = DerivativeXY(idx, idy, dxy);
00251 G4double f2xy = DerivativeXY(idx+1, idy, dxy);
00252 G4double f3xy = DerivativeXY(idx+1, idy+1, dxy);
00253 G4double f4xy = DerivativeXY(idx, idy+1, dxy);
00254
00255 cache->lastValue =
00256 f1 + f1y*h2 + (3*(f4-f1) - 2*f1y - f4y)*h22 + (2*(f1 - f4) + f1y + f4y)*h23
00257 + f1x*h1 + f1xy*h1*h2 +(3*(f4x - f1x) - 2*f1xy - f4xy)*h1*h22
00258 + (2*(f1x - f4x) + f1xy + f4xy)*h1*h23
00259 + (3*(f2 - f1) - 2*f1x - f2x)*h12 + (3*f2y - 3*f1y - 2*f1xy - f2xy)*h12*h2
00260 + (9*(f1 - f2 + f3 - f4) + 6*f1x + 3*f2x - 3*f3x - 6*f4x + 6*f1y - 6*f2y
00261 - 3*f3y + 3*f4y + 4*f1xy + 2*f2xy + f3xy + 2*f4xy)*h12*h22
00262 + (6*(-f1 + f2 - f3 + f4) - 4*f1x - 2*f2x + 2*f3x + 4*f4x - 3*f1y
00263 + 3*f2y + 3*f3y - 3*f4y - 2*f1xy - f2xy - f3xy - 2*f4xy)*h12*h23
00264 + (2*(f1 - f2) + f1x + f2x)*h13 + (2*(f1y - f2y) + f1xy + f2xy)*h13*h2
00265 + (6*(-f1 + f2 -f3 + f4) + 3*(-f1x - f2x + f3x + f4x) - 4*f1y
00266 + 4*f2y + 2*f3y - 2*f4y - 2*f1xy - 2*f2xy - f3xy - f4xy)*h13*h22
00267 + (4*(f1 - f2 + f3 - f4) + 2*(f1x + f2x - f3x - f4x)
00268 + 2*(f1y - f2y - f3y + f4y) + f1xy + f2xy + f3xy + f4xy)*h13*h23;
00269 }
00270
00271
00272
00273 void
00274 G4Physics2DVector::PutVectors(const std::vector<G4double>& vecX,
00275 const std::vector<G4double>& vecY)
00276 {
00277 ClearVectors();
00278 numberOfXNodes = vecX.size();
00279 numberOfYNodes = vecY.size();
00280 PrepareVectors();
00281 if(!cache) { cache = new G4Physics2DVectorCache(); }
00282 cache->Clear();
00283 for(size_t i = 0; i<numberOfXNodes; ++i) {
00284 xVector[i] = vecX[i];
00285 }
00286 for(size_t j = 0; j<numberOfYNodes; ++j) {
00287 yVector[j] = vecY[j];
00288 }
00289 }
00290
00291
00292
00293 void G4Physics2DVector::Store(std::ofstream& out)
00294 {
00295
00296 G4int prec = out.precision();
00297 out << G4int(type) << " " << numberOfXNodes << " " << numberOfYNodes
00298 << G4endl;
00299 out << std::setprecision(5);
00300
00301
00302 for(size_t i = 0; i<numberOfXNodes-1; ++i) {
00303 out << xVector[i] << " ";
00304 }
00305 out << xVector[numberOfXNodes-1] << G4endl;
00306 for(size_t j = 0; j<numberOfYNodes-1; ++j) {
00307 out << yVector[j] << " ";
00308 }
00309 out << yVector[numberOfYNodes-1] << G4endl;
00310 for(size_t j = 0; j<numberOfYNodes; ++j) {
00311 for(size_t i = 0; i<numberOfXNodes-1; ++i) {
00312 out << GetValue(i, j) << " ";
00313 }
00314 out << GetValue(numberOfXNodes-1,j) << G4endl;
00315 }
00316 out.precision(prec);
00317 out.close();
00318 }
00319
00320
00321
00322 G4bool G4Physics2DVector::Retrieve(std::ifstream& in)
00323 {
00324
00325 cache->Clear();
00326 ClearVectors();
00327
00328
00329 G4int k;
00330 in >> k >> numberOfXNodes >> numberOfYNodes;
00331 if (in.fail()) { return false; }
00332 PrepareVectors();
00333 type = G4PhysicsVectorType(k);
00334
00335
00336 G4double val;
00337 for(size_t i = 0; i<numberOfXNodes; ++i) {
00338 in >> xVector[i];
00339 if (in.fail()) { return false; }
00340 }
00341 for(size_t j = 0; j<numberOfYNodes; ++j) {
00342 in >> yVector[j];
00343 if (in.fail()) { return false; }
00344 }
00345 for(size_t j = 0; j<numberOfYNodes; ++j) {
00346 for(size_t i = 0; i<numberOfXNodes; ++i) {
00347 in >> val;
00348 if (in.fail()) { return false; }
00349 PutValue(i, j, val);
00350 }
00351 }
00352 in.close();
00353 return true;
00354 }
00355
00356
00357
00358 void
00359 G4Physics2DVector::ScaleVector(G4double factor)
00360 {
00361 G4double val;
00362 for(size_t j = 0; j<numberOfYNodes; ++j) {
00363 for(size_t i = 0; i<numberOfXNodes; ++i) {
00364 val = GetValue(i, j)*factor;
00365 PutValue(i, j, val);
00366 }
00367 }
00368 }
00369
00370
00371
00372 size_t
00373 G4Physics2DVector::FindBinLocation(G4double z,
00374 const G4PV2DDataVector& v)
00375 {
00376 size_t lowerBound = 0;
00377 size_t upperBound = v.size() - 2;
00378
00379 while (lowerBound <= upperBound)
00380 {
00381 size_t midBin = (lowerBound + upperBound)/2;
00382 if( z < v[midBin] ) { upperBound = midBin-1; }
00383 else { lowerBound = midBin+1; }
00384 }
00385
00386 return upperBound;
00387 }
00388
00389