Geant4  10.00.p03
G4PhysicsVector.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 // $Id: G4PhysicsVector.cc 74256 2013-10-02 14:24:02Z gcosmo $
28 //
29 //
30 // --------------------------------------------------------------
31 // GEANT 4 class implementation file
32 //
33 // G4PhysicsVector.cc
34 //
35 // History:
36 // 02 Dec. 1995, G.Cosmo : Structure created based on object model
37 // 03 Mar. 1996, K.Amako : Implemented the 1st version
38 // 01 Jul. 1996, K.Amako : Hidden bin from the user introduced
39 // 12 Nov. 1998, K.Amako : A bug in GetVectorLength() fixed
40 // 11 Nov. 2000, H.Kurashige : use STL vector for dataVector and binVector
41 // 18 Jan. 2001, H.Kurashige : removed ptrNextTable
42 // 09 Mar. 2001, H.Kurashige : added G4PhysicsVector type
43 // 05 Sep. 2008, V.Ivanchenko : added protections for zero-length vector
44 // 11 May 2009, A.Bagulya : added new implementation of methods
45 // ComputeSecondDerivatives - first derivatives at edge points
46 // should be provided by a user
47 // FillSecondDerivatives - default computation base on "not-a-knot"
48 // algorithm
49 // 19 Jun. 2009, V.Ivanchenko : removed hidden bin
50 // 17 Nov. 2009, H.Kurashige : use pointer for DataVector
51 // 04 May 2010 H.Kurashige : use G4PhyscisVectorCache
52 // 28 May 2010 H.Kurashige : Stop using pointers to G4PVDataVector
53 // 16 Aug. 2011 H.Kurashige : Add dBin, baseBin and verboseLevel
54 // --------------------------------------------------------------
55 
56 #include <iomanip>
57 #include "G4PhysicsVector.hh"
58 
60 
61 // --------------------------------------------------------------
62 
64  : type(T_G4PhysicsVector),
65  edgeMin(0.), edgeMax(0.), numberOfNodes(0),
66  useSpline(false),
67  dBin(0.), baseBin(0.),
68  verboseLevel(0)
69 {
71  // g4pow = G4Pow::GetInstance();
72 }
73 
74 // --------------------------------------------------------------
75 
77 {
78 }
79 
80 // --------------------------------------------------------------
81 
83 {
84  // g4pow = G4Pow::GetInstance();
85 
86  dBin = right.dBin;
87  baseBin = right.baseBin;
88  verboseLevel = right.verboseLevel;
89 
90  DeleteData();
91  CopyData(right);
92 }
93 
94 // --------------------------------------------------------------
95 
97 {
98  if (&right==this) { return *this; }
99  dBin = right.dBin;
100  baseBin = right.baseBin;
101  verboseLevel = right.verboseLevel;
102 
103  DeleteData();
104  CopyData(right);
105  return *this;
106 }
107 
108 // --------------------------------------------------------------
109 
111 {
112  return (this == &right);
113 }
114 
115 // --------------------------------------------------------------
116 
118 {
119  return (this != &right);
120 }
121 
122 // --------------------------------------------------------------
123 
125 {
126  useSpline = false;
127  secDerivative.clear();
128 }
129 
130 // --------------------------------------------------------------
131 
133 {
134  type = vec.type;
135  edgeMin = vec.edgeMin;
136  edgeMax = vec.edgeMax;
138  useSpline = vec.useSpline;
139 
140  size_t i;
141  dataVector.clear();
142  for(i=0; i<(vec.dataVector).size(); i++){
143  dataVector.push_back( (vec.dataVector)[i] );
144  }
145  binVector.clear();
146  for(i=0; i<(vec.binVector).size(); i++){
147  binVector.push_back( (vec.binVector)[i] );
148  }
149  secDerivative.clear();
150  for(i=0; i<(vec.secDerivative).size(); i++){
151  secDerivative.push_back( (vec.secDerivative)[i] );
152  }
153 }
154 
155 // --------------------------------------------------------------
156 
158 {
159  return binVector[binNumber];
160 }
161 
162 // --------------------------------------------------------------
163 
164 G4bool G4PhysicsVector::Store(std::ofstream& fOut, G4bool ascii)
165 {
166  // Ascii mode
167  if (ascii)
168  {
169  fOut << *this;
170  return true;
171  }
172  // Binary Mode
173 
174  // binning
175  fOut.write((char*)(&edgeMin), sizeof edgeMin);
176  fOut.write((char*)(&edgeMax), sizeof edgeMax);
177  fOut.write((char*)(&numberOfNodes), sizeof numberOfNodes);
178 
179  // contents
180  size_t size = dataVector.size();
181  fOut.write((char*)(&size), sizeof size);
182 
183  G4double* value = new G4double[2*size];
184  for(size_t i = 0; i < size; ++i)
185  {
186  value[2*i] = binVector[i];
187  value[2*i+1]= dataVector[i];
188  }
189  fOut.write((char*)(value), 2*size*(sizeof (G4double)));
190  delete [] value;
191 
192  return true;
193 }
194 
195 // --------------------------------------------------------------
196 
197 G4bool G4PhysicsVector::Retrieve(std::ifstream& fIn, G4bool ascii)
198 {
199  // clear properties;
200  dataVector.clear();
201  binVector.clear();
202  secDerivative.clear();
203 
204  // retrieve in ascii mode
205  if (ascii){
206  // binning
207  fIn >> edgeMin >> edgeMax >> numberOfNodes;
208  if (fIn.fail()) { return false; }
209  // contents
210  G4int siz=0;
211  fIn >> siz;
212  if (fIn.fail()) { return false; }
213  if (siz<=0)
214  {
215 #ifdef G4VERBOSE
216  G4cerr << "G4PhysicsVector::Retrieve():";
217  G4cerr << " Invalid vector size: " << siz << G4endl;
218 #endif
219  return false;
220  }
221 
222  binVector.reserve(siz);
223  dataVector.reserve(siz);
224  G4double vBin, vData;
225 
226  for(G4int i = 0; i < siz ; i++)
227  {
228  vBin = 0.;
229  vData= 0.;
230  fIn >> vBin >> vData;
231  if (fIn.fail()) { return false; }
232  binVector.push_back(vBin);
233  dataVector.push_back(vData);
234  }
235 
236  // to remove any inconsistency
237  numberOfNodes = siz;
238  edgeMin = binVector[0];
239  edgeMax = binVector[numberOfNodes-1];
240  return true ;
241  }
242 
243  // retrieve in binary mode
244  // binning
245  fIn.read((char*)(&edgeMin), sizeof edgeMin);
246  fIn.read((char*)(&edgeMax), sizeof edgeMax);
247  fIn.read((char*)(&numberOfNodes), sizeof numberOfNodes );
248 
249  // contents
250  size_t size;
251  fIn.read((char*)(&size), sizeof size);
252 
253  G4double* value = new G4double[2*size];
254  fIn.read((char*)(value), 2*size*(sizeof(G4double)) );
255  if (G4int(fIn.gcount()) != G4int(2*size*(sizeof(G4double))) )
256  {
257  delete [] value;
258  return false;
259  }
260 
261  binVector.reserve(size);
262  dataVector.reserve(size);
263  for(size_t i = 0; i < size; ++i)
264  {
265  binVector.push_back(value[2*i]);
266  dataVector.push_back(value[2*i+1]);
267  }
268  delete [] value;
269 
270  // to remove any inconsistency
271  numberOfNodes = size;
272  edgeMin = binVector[0];
274 
275  return true;
276 }
277 
278 // --------------------------------------------------------------
279 
280 void
282 {
283  size_t n = dataVector.size();
284  size_t i;
285  if(n > 0) {
286  for(i=0; i<n; ++i) {
287  binVector[i] *= factorE;
288  dataVector[i] *= factorV;
289  }
290  }
291  // n = secDerivative.size();
292  // if(n > 0) { for(i=0; i<n; ++i) { secDerivative[i] *= factorV; } }
293  secDerivative.clear();
294 
295  edgeMin *= factorE;
296  edgeMax *= factorE;
297 }
298 
299 // --------------------------------------------------------------
300 
301 void
303  G4double endPointDerivative)
304  // A standard method of computation of second derivatives
305  // First derivatives at the first and the last point should be provided
306  // See for example W.H. Press et al. "Numerical recipes in C"
307  // Cambridge University Press, 1997.
308 {
309  if(4 > numberOfNodes) // cannot compute derivatives for less than 4 bins
310  {
312  return;
313  }
314 
315  if(!SplinePossible()) { return; }
316 
317  useSpline = true;
318 
319  G4int n = numberOfNodes-1;
320 
321  G4double* u = new G4double [n];
322 
323  G4double p, sig, un;
324 
325  u[0] = (6.0/(binVector[1]-binVector[0]))
326  * ((dataVector[1]-dataVector[0])/(binVector[1]-binVector[0])
327  - firstPointDerivative);
328 
329  secDerivative[0] = - 0.5;
330 
331  // Decomposition loop for tridiagonal algorithm. secDerivative[i]
332  // and u[i] are used for temporary storage of the decomposed factors.
333 
334  for(G4int i=1; i<n; ++i)
335  {
336  sig = (binVector[i]-binVector[i-1]) / (binVector[i+1]-binVector[i-1]);
337  p = sig*(secDerivative[i-1]) + 2.0;
338  secDerivative[i] = (sig - 1.0)/p;
339  u[i] = (dataVector[i+1]-dataVector[i])/(binVector[i+1]-binVector[i])
340  - (dataVector[i]-dataVector[i-1])/(binVector[i]-binVector[i-1]);
341  u[i] = 6.0*u[i]/(binVector[i+1]-binVector[i-1]) - sig*u[i-1]/p;
342  }
343 
344  sig = (binVector[n-1]-binVector[n-2]) / (binVector[n]-binVector[n-2]);
345  p = sig*secDerivative[n-2] + 2.0;
346  un = (6.0/(binVector[n]-binVector[n-1]))
347  *(endPointDerivative -
348  (dataVector[n]-dataVector[n-1])/(binVector[n]-binVector[n-1])) - u[n-1]/p;
349  secDerivative[n] = un/(secDerivative[n-1] + 2.0);
350 
351  // The back-substitution loop for the triagonal algorithm of solving
352  // a linear system of equations.
353 
354  for(G4int k=n-1; k>0; --k)
355  {
356  secDerivative[k] *=
357  (secDerivative[k+1] -
358  u[k]*(binVector[k+1]-binVector[k-1])/(binVector[k+1]-binVector[k]));
359  }
360  secDerivative[0] = 0.5*(u[0] - secDerivative[1]);
361 
362  delete [] u;
363 }
364 
365 // --------------------------------------------------------------
366 
368  // Computation of second derivatives using "Not-a-knot" endpoint conditions
369  // B.I. Kvasov "Methods of shape-preserving spline approximation"
370  // World Scientific, 2000
371 {
372  if(5 > numberOfNodes) // cannot compute derivatives for less than 4 points
373  {
375  return;
376  }
377 
378  if(!SplinePossible()) { return; }
379 
380  useSpline = true;
381 
382  G4int n = numberOfNodes-1;
383 
384  G4double* u = new G4double [n];
385 
386  G4double p, sig;
387 
388  u[1] = ((dataVector[2]-dataVector[1])/(binVector[2]-binVector[1]) -
389  (dataVector[1]-dataVector[0])/(binVector[1]-binVector[0]));
390  u[1] = 6.0*u[1]*(binVector[2]-binVector[1])
391  / ((binVector[2]-binVector[0])*(binVector[2]-binVector[0]));
392 
393  // Decomposition loop for tridiagonal algorithm. secDerivative[i]
394  // and u[i] are used for temporary storage of the decomposed factors.
395 
396  secDerivative[1] = (2.0*binVector[1]-binVector[0]-binVector[2])
397  / (2.0*binVector[2]-binVector[0]-binVector[1]);
398 
399  for(G4int i=2; i<n-1; ++i)
400  {
401  sig = (binVector[i]-binVector[i-1]) / (binVector[i+1]-binVector[i-1]);
402  p = sig*secDerivative[i-1] + 2.0;
403  secDerivative[i] = (sig - 1.0)/p;
404  u[i] = (dataVector[i+1]-dataVector[i])/(binVector[i+1]-binVector[i])
405  - (dataVector[i]-dataVector[i-1])/(binVector[i]-binVector[i-1]);
406  u[i] = (6.0*u[i]/(binVector[i+1]-binVector[i-1])) - sig*u[i-1]/p;
407  }
408 
409  sig = (binVector[n-1]-binVector[n-2]) / (binVector[n]-binVector[n-2]);
410  p = sig*secDerivative[n-3] + 2.0;
411  u[n-1] = (dataVector[n]-dataVector[n-1])/(binVector[n]-binVector[n-1])
412  - (dataVector[n-1]-dataVector[n-2])/(binVector[n-1]-binVector[n-2]);
413  u[n-1] = 6.0*sig*u[n-1]/(binVector[n]-binVector[n-2])
414  - (2.0*sig - 1.0)*u[n-2]/p;
415 
416  p = (1.0+sig) + (2.0*sig-1.0)*secDerivative[n-2];
417  secDerivative[n-1] = u[n-1]/p;
418 
419  // The back-substitution loop for the triagonal algorithm of solving
420  // a linear system of equations.
421 
422  for(G4int k=n-2; k>1; --k)
423  {
424  secDerivative[k] *=
425  (secDerivative[k+1] -
426  u[k]*(binVector[k+1]-binVector[k-1])/(binVector[k+1]-binVector[k]));
427  }
428  secDerivative[n] = (secDerivative[n-1] - (1.0-sig)*secDerivative[n-2])/sig;
429  sig = 1.0 - ((binVector[2]-binVector[1])/(binVector[2]-binVector[0]));
430  secDerivative[1] *= (secDerivative[2] - u[1]/(1.0-sig));
431  secDerivative[0] = (secDerivative[1] - sig*secDerivative[2])/(1.0-sig);
432 
433  delete [] u;
434 }
435 
436 // --------------------------------------------------------------
437 
438 void
440  // A simplified method of computation of second derivatives
441 {
442  if(3 > numberOfNodes) // cannot compute derivatives for less than 4 bins
443  {
444  useSpline = false;
445  return;
446  }
447 
448  if(!SplinePossible()) { return; }
449 
450  useSpline = true;
451 
452  size_t n = numberOfNodes-1;
453 
454  for(size_t i=1; i<n; ++i)
455  {
456  secDerivative[i] =
457  3.0*((dataVector[i+1]-dataVector[i])/(binVector[i+1]-binVector[i]) -
458  (dataVector[i]-dataVector[i-1])/(binVector[i]-binVector[i-1]))
459  /(binVector[i+1]-binVector[i-1]);
460  }
461  secDerivative[n] = secDerivative[n-1];
463 }
464 
465 // --------------------------------------------------------------
466 
468  // Initialise second derivative array. If neighbor energy coincide
469  // or not ordered than spline cannot be applied
470 {
471  G4bool result = true;
472  secDerivative.clear();
473  secDerivative.reserve(numberOfNodes);
474  for(size_t j=0; j<numberOfNodes; ++j)
475  {
476  secDerivative.push_back(0.0);
477  if(j > 0)
478  {
479  if(binVector[j]-binVector[j-1] <= 0.) {
480  result = false;
481  secDerivative.clear();
482  break;
483  }
484  }
485  }
486  return result;
487 }
488 
489 // --------------------------------------------------------------
490 
491 std::ostream& operator<<(std::ostream& out, const G4PhysicsVector& pv)
492 {
493  // binning
494  out << std::setprecision(12) << pv.edgeMin << " "
495  << pv.edgeMax << " " << pv.numberOfNodes << G4endl;
496 
497  // contents
498  out << pv.dataVector.size() << G4endl;
499  for(size_t i = 0; i < pv.dataVector.size(); i++)
500  {
501  out << pv.binVector[i] << " " << pv.dataVector[i] << G4endl;
502  }
503  out << std::setprecision(6);
504 
505  return out;
506 }
507 
508 //---------------------------------------------------------------
509 
510 G4double
511 G4PhysicsVector::Value(G4double theEnergy, size_t& lastIdx) const
512 {
513  G4double y;
514  if(theEnergy <= edgeMin) {
515  lastIdx = 0;
516  y = dataVector[0];
517  } else if(theEnergy >= edgeMax) {
518  lastIdx = numberOfNodes-1;
519  y = dataVector[lastIdx];
520  } else {
521  lastIdx = FindBin(theEnergy, lastIdx);
522  y = Interpolation(lastIdx, theEnergy);
523  }
524  return y;
525 }
526 
527 //---------------------------------------------------------------
528 
530 {
531  if(1 >= numberOfNodes) { return 0.0; }
532  size_t n1 = 0;
533  size_t n2 = numberOfNodes/2;
534  size_t n3 = numberOfNodes - 1;
535  G4double y = rand*dataVector[n3];
536  while (n1 + 1 != n3)
537  {
538  if (y > dataVector[n2])
539  { n1 = n2; }
540  else
541  { n3 = n2; }
542  n2 = (n3 + n1 + 1)/2;
543  }
544  G4double res = binVector[n1];
545  G4double del = dataVector[n3] - dataVector[n1];
546  if(del > 0.0) {
547  res += (y - dataVector[n1])*(binVector[n3] - res)/del;
548  }
549  return res;
550 }
551 
552 //---------------------------------------------------------------
G4PVDataVector dataVector
G4double FindLinearEnergy(G4double rand) const
G4double Interpolation(size_t idx, G4double energy) const
void CopyData(const G4PhysicsVector &vec)
G4ThreadLocal G4Allocator< G4PhysicsVector > * fpPVAllocator
G4PhysicsVector(G4bool spline=false)
void ComputeSecondDerivatives(G4double firstPointDerivative, G4double endPointDerivative)
G4double GetLowEdgeEnergy(size_t binNumber) const
#define G4ThreadLocal
Definition: tls.hh:52
int G4int
Definition: G4Types.hh:78
void FillSecondDerivatives()
G4PVDataVector binVector
G4int operator!=(const G4PhysicsVector &right) const
std::ostream & operator<<(std::ostream &out, const G4PhysicsVector &pv)
G4int operator==(const G4PhysicsVector &right) const
bool G4bool
Definition: G4Types.hh:79
G4PhysicsVectorType type
G4PVDataVector secDerivative
virtual void ScaleVector(G4double factorE, G4double factorV)
size_t FindBin(G4double energy, size_t idx) const
const G4int n
G4double Value(G4double theEnergy, size_t &lastidx) const
virtual G4bool Store(std::ofstream &fOut, G4bool ascii=false)
virtual G4bool Retrieve(std::ifstream &fIn, G4bool ascii=false)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void ComputeSecDerivatives()
virtual ~G4PhysicsVector()
G4PhysicsVector & operator=(const G4PhysicsVector &)
G4GLOB_DLL std::ostream G4cerr