Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/global/management/src/G4Physics2DVector.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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