G4Physics2DVector.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id:$
00028 // GEANT4 tag $Name: not supported by cvs2svn $
00029 //
00030 // 
00031 // --------------------------------------------------------------
00032 //      GEANT 4 class implementation file
00033 //
00034 //  G4Physics2DVector.cc
00035 //
00036 //  Author:        Vladimir Ivanchenko
00037 //
00038 //  Creation date: 25.09.2011
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     // Bicubic interpolation according to 
00209     // 1. H.M. Antia, "Numerical Methods for Scientists and Engineers",
00210     //    MGH, 1991. 
00211     // 2. W.H. Press et al., "Numerical recipes. The Art of Scientific 
00212     //    Computing", Cambridge University Press, 2007. 
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     // Three derivatives at each of four points (1-4) defining the 
00236     // subregion are computed by numerical centered differencing from 
00237     // the functional values already tabulated on the grid. 
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   // binning
00296   G4int prec = out.precision();
00297   out << G4int(type) << " " << numberOfXNodes << " " << numberOfYNodes 
00298       << G4endl; 
00299   out << std::setprecision(5);
00300 
00301   // contents
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   // initialisation
00325   cache->Clear();
00326   ClearVectors();
00327 
00328   // binning
00329   G4int k;
00330   in >> k >> numberOfXNodes >> numberOfYNodes;
00331   if (in.fail())  { return false; }
00332   PrepareVectors();
00333   type = G4PhysicsVectorType(k); 
00334 
00335   // contents
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 // --------------------------------------------------------------

Generated on Mon May 27 17:49:19 2013 for Geant4 by  doxygen 1.4.7