Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Physics2DVector.cc
Go to the documentation of this file.
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 // --------------------------------------------------------------
27 // GEANT 4 class implementation file
28 //
29 // G4Physics2DVector.cc
30 //
31 // Author: Vladimir Ivanchenko
32 //
33 // Creation date: 25.09.2011
34 //
35 // --------------------------------------------------------------
36 
37 #include <iomanip>
38 
39 #include "G4Physics2DVector.hh"
40 
41 // --------------------------------------------------------------
42 
44  : type(T_G4PhysicsFreeVector),
45  numberOfXNodes(0), numberOfYNodes(0),
46  verboseLevel(0), useBicubic(false)
47 {}
48 
49 // --------------------------------------------------------------
50 
52  : type(T_G4PhysicsFreeVector),
53  numberOfXNodes(nx), numberOfYNodes(ny),
54  verboseLevel(0), useBicubic(false)
55 {
57 }
58 
59 // --------------------------------------------------------------
60 
62 {
63  ClearVectors();
64 }
65 
66 // --------------------------------------------------------------
67 
69 {
70  type = right.type;
71 
72  numberOfXNodes = right.numberOfXNodes;
73  numberOfYNodes = right.numberOfYNodes;
74 
75  verboseLevel = right.verboseLevel;
76  useBicubic = right.useBicubic;
77 
78  xVector = right.xVector;
79  yVector = right.yVector;
80 
82  CopyData(right);
83 }
84 
85 // --------------------------------------------------------------
86 
88 {
89  if (&right==this) { return *this; }
90  ClearVectors();
91 
92  type = right.type;
93 
94  numberOfXNodes = right.numberOfXNodes;
95  numberOfYNodes = right.numberOfYNodes;
96 
97  verboseLevel = right.verboseLevel;
98  useBicubic = right.useBicubic;
99 
100  PrepareVectors();
101  CopyData(right);
102 
103  return *this;
104 }
105 
106 // --------------------------------------------------------------
107 
109 {
110  xVector.resize(numberOfXNodes,0.);
111  yVector.resize(numberOfYNodes,0.);
112  value.resize(numberOfYNodes,0);
113  for(size_t j=0; j<numberOfYNodes; ++j) {
115  v->resize(numberOfXNodes,0.);
116  value[j] = v;
117  }
118 }
119 
120 // --------------------------------------------------------------
121 
123 {
124  for(size_t j=0; j<numberOfYNodes; ++j) {
125  delete value[j];
126  }
127 }
128 
129 // --------------------------------------------------------------
130 
132 {
133  for(size_t i=0; i<numberOfXNodes; ++i) {
134  xVector[i] = right.xVector[i];
135  }
136  for(size_t j=0; j<numberOfYNodes; ++j) {
137  yVector[j] = right.yVector[j];
138  G4PV2DDataVector* v0 = right.value[j];
139  for(size_t i=0; i<numberOfXNodes; ++i) {
140  PutValue(i,j,(*v0)[i]);
141  }
142  }
143 }
144 
145 // --------------------------------------------------------------
146 
148  size_t& idx, size_t& idy) const
149 {
150  G4double x = xx;
151  G4double y = yy;
152 
153  // no interpolation outside the table
154  if(x < xVector[0]) {
155  x = xVector[0];
156  } else if(x > xVector[numberOfXNodes - 1]) {
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  }
164 
165  // find bins
166  idx = FindBinLocationX(x, idx);
167  idy = FindBinLocationY(y, idy);
168 
169  // interpolate
170  if(useBicubic) {
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));
183  }
184 }
185 
186 // --------------------------------------------------------------
187 
188 G4double
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)
258 {
259  ClearVectors();
260  numberOfXNodes = vecX.size();
261  numberOfYNodes = vecY.size();
262  PrepareVectors();
263  for(size_t i = 0; i<numberOfXNodes; ++i) {
264  xVector[i] = vecX[i];
265  }
266  for(size_t j = 0; j<numberOfYNodes; ++j) {
267  yVector[j] = vecY[j];
268  }
269 }
270 
271 // --------------------------------------------------------------
272 
273 void G4Physics2DVector::Store(std::ofstream& out) const
274 {
275  // binning
276  G4int prec = out.precision();
277  out << G4int(type) << " " << numberOfXNodes << " " << numberOfYNodes
278  << G4endl;
279  out << std::setprecision(5);
280 
281  // contents
282  for(size_t i = 0; i<numberOfXNodes-1; ++i) {
283  out << xVector[i] << " ";
284  }
285  out << xVector[numberOfXNodes-1] << G4endl;
286  for(size_t j = 0; j<numberOfYNodes-1; ++j) {
287  out << yVector[j] << " ";
288  }
289  out << yVector[numberOfYNodes-1] << G4endl;
290  for(size_t j = 0; j<numberOfYNodes; ++j) {
291  for(size_t i = 0; i<numberOfXNodes-1; ++i) {
292  out << GetValue(i, j) << " ";
293  }
294  out << GetValue(numberOfXNodes-1,j) << G4endl;
295  }
296  out.precision(prec);
297  out.close();
298 }
299 
300 // --------------------------------------------------------------
301 
303 {
304  // initialisation
305  ClearVectors();
306 
307  // binning
308  G4int k;
309  in >> k >> numberOfXNodes >> numberOfYNodes;
310  if (in.fail() || 0 >= numberOfXNodes || 0 >= numberOfYNodes ||
311  numberOfXNodes >= INT_MAX || numberOfYNodes >= INT_MAX) {
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;
319  }
320  PrepareVectors();
321  type = G4PhysicsVectorType(k);
322 
323  // contents
324  G4double val;
325  for(size_t i = 0; i<numberOfXNodes; ++i) {
326  in >> xVector[i];
327  if (in.fail()) { return false; }
328  }
329  for(size_t j = 0; j<numberOfYNodes; ++j) {
330  in >> yVector[j];
331  if (in.fail()) { return false; }
332  }
333  for(size_t j = 0; j<numberOfYNodes; ++j) {
334  for(size_t i = 0; i<numberOfXNodes; ++i) {
335  in >> val;
336  if (in.fail()) { return false; }
337  PutValue(i, j, val);
338  }
339  }
340  in.close();
341  return true;
342 }
343 
344 // --------------------------------------------------------------
345 
346 void
348 {
349  G4double val;
350  for(size_t j = 0; j<numberOfYNodes; ++j) {
351  for(size_t i = 0; i<numberOfXNodes; ++i) {
352  val = GetValue(i, j)*factor;
353  PutValue(i, j, val);
354  }
355  }
356 }
357 
358 // --------------------------------------------------------------
359 
360 size_t
362  const G4PV2DDataVector& v) const
363 {
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 
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  }
388 
389  // find bins
390  idy = FindBinLocationY(y, idy);
391 
392  G4double x1 = InterpolateLinearX(*(value[idy]), rand);
393  G4double x2 = InterpolateLinearX(*(value[idy+1]), rand);
394  G4double res = x1;
395  G4double del = yVector[idy+1] - yVector[idy];
396  if(del != 0.0) {
397  res += (x2 - x1)*(y - yVector[idy])/del;
398  }
399  return res;
400 }
401 
402 // --------------------------------------------------------------
403 
404 G4double G4Physics2DVector::InterpolateLinearX(G4PV2DDataVector& v,
405  G4double rand) const
406 {
407  size_t nn = v.size();
408  if(1 >= nn) { return 0.0; }
409  size_t n1 = 0;
410  size_t n2 = nn/2;
411  size_t n3 = nn - 1;
412  G4double y = rand*v[n3];
413  while (n1 + 1 != n3)
414  {
415  if (y > v[n2])
416  { n1 = n2; }
417  else
418  { n3 = n2; }
419  n2 = (n3 + n1 + 1)/2;
420  }
421  G4double res = xVector[n1];
422  G4double del = v[n3] - v[n1];
423  if(del > 0.0) {
424  res += (y - v[n1])*(xVector[n3] - res)/del;
425  }
426  return res;
427 }
428 
429 // --------------------------------------------------------------
std::vector< G4double > G4PV2DDataVector
G4double GetValue(size_t idx, size_t idy) const
tuple bin
Definition: plottest35.py:22
G4PhysicsVectorType
size_t FindBinLocation(G4double z, const G4PV2DDataVector &) const
tuple x
Definition: test.py:50
size_t FindBinLocationX(G4double x, size_t lastidx) const
int G4int
Definition: G4Types.hh:78
void CopyData(const G4Physics2DVector &vec)
void PutValue(size_t idx, size_t idy, G4double value)
static const double prec
Definition: RanecuEngine.cc:58
G4double Value(G4double x, G4double y, size_t &lastidx, size_t &lastidy) const
const XML_Char int const XML_Char * value
Definition: expat.h:331
bool G4bool
Definition: G4Types.hh:79
G4Physics2DVector & operator=(const G4Physics2DVector &)
tuple v
Definition: test.py:18
void Store(std::ofstream &fOut) const
#define INT_MAX
Definition: templates.hh:111
G4bool Retrieve(std::ifstream &fIn)
tuple z
Definition: test.py:28
G4double BicubicInterpolation(G4double x, G4double y, size_t idx, size_t idy) const
size_t FindBinLocationY(G4double y, size_t lastidy) const
#define G4endl
Definition: G4ios.hh:61
void PutVectors(const std::vector< G4double > &vecX, const std::vector< G4double > &vecY)
double G4double
Definition: G4Types.hh:76
void ScaleVector(G4double factor)
G4double FindLinearX(G4double rand, G4double y, size_t &lastidy) const