Geant4  10.01.p02
G4PhysicsVector.icc
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.icc 83009 2014-07-24 14:51:29Z gcosmo $
28 //
29 //
30 //---------------------------------------------------------------
31 // GEANT 4 class source file
32 //
33 // G4PhysicsVector.icc
34 //
35 // Description:
36 // A physics vector which has values of energy-loss, cross-section,
37 // and other physics values of a particle in matter in a given
38 // range of the energy, momentum, etc.
39 // This class serves as the base class for a vector having various
40 // energy scale, for example like 'log', 'linear', 'free', etc.
41 //
42 //---------------------------------------------------------------
43 
44 //extern G4GLOB_DLL G4ThreadLocal G4Allocator<G4PhysicsVector> *fpPVAllocator;
45 
46 //---------------------------------------------------------------
47 /*
48 inline
49  void* G4PhysicsVector::operator new(size_t)
50 {
51  if (!fpPVAllocator) fpPVAllocator = new G4Allocator<G4PhysicsVector>;
52  return (void*)fpPVAllocator->MallocSingle();
53 }
54 
55 //---------------------------------------------------------------
56 
57 inline
58  void G4PhysicsVector::operator delete(void* aVector)
59 {
60  fpPVAllocator->FreeSingle((G4PhysicsVector*)aVector);
61 }
62 */
63 //---------------------------------------------------------------
64 
65 inline
66  G4double G4PhysicsVector::operator[](const size_t binNumber) const
67 {
68  return dataVector[binNumber];
69 }
70 
71 //---------------------------------------------------------------
72 
73 inline
74  G4double G4PhysicsVector::operator()(const size_t binNumber) const
75 {
76  return dataVector[binNumber];
77 }
78 
79 //---------------------------------------------------------------
80 
81 inline
82  G4double G4PhysicsVector::Energy(const size_t binNumber) const
83 {
84  return binVector[binNumber];
85 }
86 
87 //---------------------------------------------------------------
88 
89 inline
90  G4double G4PhysicsVector::GetMaxEnergy() const
91 {
92  return edgeMax;
93 }
94 
95 //---------------------------------------------------------------
96 
97 inline
98  size_t G4PhysicsVector::GetVectorLength() const
99 {
100  return numberOfNodes;
101 }
102 
103 //---------------------------------------------------------------
104 
105 inline
106  G4double G4PhysicsVector::GetValue(G4double theEnergy, G4bool&) const
107 {
108  size_t idx=0;
109  return Value(theEnergy, idx);
110 }
111 
112 //------------------------------------------------
113 
114 inline
115  G4double G4PhysicsVector::LinearInterpolation(size_t idx, G4double e) const
116 {
117  // Linear interpolation is used to get the value. Before this method
118  // is called it is ensured that the energy is inside the bin
119  // 0 < idx < numberOfNodes-1
120 
121  return dataVector[idx] +
122  ( dataVector[idx + 1]-dataVector[idx] ) * (e - binVector[idx])
123  /( binVector[idx + 1]-binVector[idx] );
124 }
125 
126 //---------------------------------------------------------------
127 
128 inline
129  G4double G4PhysicsVector::SplineInterpolation(size_t idx, G4double e) const
130 {
131  // Spline interpolation is used to get the value. Before this method
132  // is called it is ensured that the energy is inside the bin
133  // 0 < idx < numberOfNodes-1
134 
135  static const G4double onesixth = 1.0/6.0;
136 
137  // check bin value
138  G4double x1 = binVector[idx];
139  G4double x2 = binVector[idx + 1];
140  G4double delta = x2 - x1;
141 
142  G4double a = (x2 - e)/delta;
143  G4double b = (e - x1)/delta;
144 
145  // Final evaluation of cubic spline polynomial for return
146  G4double y1 = dataVector[idx];
147  G4double y2 = dataVector[idx + 1];
148 
149  G4double res = a*y1 + b*y2 +
150  ( (a*a*a - a)*secDerivative[idx] +
151  (b*b*b - b)*secDerivative[idx + 1] )*delta*delta*onesixth;
152 
153  return res;
154 }
155 
156 //---------------------------------------------------------------
157 
158 inline
159  G4double G4PhysicsVector::Interpolation(size_t idx, G4double e) const
160 {
161  G4double res;
162  if(useSpline) { res = SplineInterpolation(idx, e); }
163  else { res = LinearInterpolation(idx, e); }
164  return res;
165 }
166 
167 //---------------------------------------------------------------
168 
169 inline
170  void G4PhysicsVector::PutValue(size_t binNumber, G4double theValue)
171 {
172  dataVector[binNumber] = theValue;
173 }
174 
175 //---------------------------------------------------------------
176 
177 inline
178  G4bool G4PhysicsVector::IsFilledVectorExist() const
179 {
180  G4bool status=false;
181 
182  if(numberOfNodes > 0) { status=true; }
183  return status;
184 }
185 
186 //---------------------------------------------------------------
187 
188 inline
189  G4PhysicsVectorType G4PhysicsVector::GetType() const
190 {
191  return type;
192 }
193 
194 //---------------------------------------------------------------
195 
196 // Flag useSpline is "true" only if second derivatives are filled
197 inline
198  void G4PhysicsVector::SetSpline(G4bool val)
199 {
200  if(val) {
201  if(0 == secDerivative.size() && 0 < dataVector.size()) {
202  FillSecondDerivatives();
203  }
204  } else {
205  useSpline = false;
206  secDerivative.clear();
207  }
208 }
209 
210 //---------------------------------------------------------------
211 
212 inline
213 void G4PhysicsVector::SetVerboseLevel(G4int value)
214 {
215  verboseLevel = value;
216 }
217 
218 //---------------------------------------------------------------
219 
220 inline
221 G4int G4PhysicsVector::GetVerboseLevel(G4int)
222 {
223  return verboseLevel;
224 }
225 
226 //---------------------------------------------------------------
227 
228 inline
229 size_t G4PhysicsVector::FindBinLocation(G4double theEnergy) const
230 {
231  size_t bin;
232  if(type == T_G4PhysicsLogVector) {
233  bin = size_t(G4Log(theEnergy)/dBin - baseBin);
234  if(bin + 2 > numberOfNodes) { bin = numberOfNodes - 2; }
235  else if(bin > 0 && theEnergy < binVector[bin]) { --bin; }
236  else if(bin + 2 < numberOfNodes && theEnergy > binVector[bin+1])
237  { ++bin; }
238  } else if(type == T_G4PhysicsLinearVector) {
239  bin = size_t( theEnergy/dBin - baseBin );
240  if(bin + 2 > numberOfNodes) { bin = numberOfNodes - 2; }
241  else if(bin > 0 && theEnergy < binVector[bin]) { --bin; }
242  else if(bin + 2 < numberOfNodes && theEnergy > binVector[bin+1])
243  { ++bin; }
244  } else {
245  bin = 0;
246  size_t bin2;
247  size_t bin3 = numberOfNodes - 1;
248  while (bin != bin3 - 1) {
249  bin2 = bin + (bin3 - bin + 1)/2;
250  if (theEnergy > binVector[bin2]) { bin = bin2; }
251  else { bin3 = bin2; }
252  }
253 /*
254 // Bin location proposed by K.Genser (FNAL)
255 // V.I. Usage of this algorithm provides identical results
256 // no CPU advantage is observed in EM tests
257 // if new validation information will be known this code may be used
258  G4PVDataVector::const_iterator it =
259  std::lower_bound(binVector.begin(), binVector.end(), theEnergy);
260  bin = it - binVector.begin() - 1;
261 */
262  }
263  return bin;
264 }
265 
266 //---------------------------------------------------------------
267 
268 inline size_t G4PhysicsVector::FindBin(G4double e, size_t idx) const
269 {
270  size_t id = idx;
271  if(e < binVector[1]) {
272  id = 0;
273  } else if(e >= binVector[numberOfNodes-2]) {
274  id = numberOfNodes - 2;
275  } else if(idx >= numberOfNodes || e < binVector[idx]
276  || e > binVector[idx+1]) {
277  id = FindBinLocation(e);
278  }
279  return id;
280 }
281 
282 //---------------------------------------------------------------
283 
284 inline
285  G4double G4PhysicsVector::Value(G4double theEnergy) const
286 {
287  size_t idx=0;
288  return Value(theEnergy, idx);
289 }
290 
291 //---------------------------------------------------------------