Geant4  10.00.p01
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 79149 2014-02-19 15:08:06Z 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::Value(G4double theEnergy) const
67 {
68  size_t idx=0;
69  return Value(theEnergy, idx);
70 }
71 
72 //---------------------------------------------------------------
73 
74 inline
75  G4double G4PhysicsVector::operator[](const size_t binNumber) const
76 {
77  return dataVector[binNumber];
78 }
79 
80 //---------------------------------------------------------------
81 
82 inline
83  G4double G4PhysicsVector::operator()(const size_t binNumber) const
84 {
85  return dataVector[binNumber];
86 }
87 
88 //---------------------------------------------------------------
89 
90 inline
91  G4double G4PhysicsVector::Energy(const size_t binNumber) const
92 {
93  return binVector[binNumber];
94 }
95 
96 //---------------------------------------------------------------
97 
98 inline
99  G4double G4PhysicsVector::GetMaxEnergy() const
100 {
101  return edgeMax;
102 }
103 
104 //---------------------------------------------------------------
105 
106 inline
107  size_t G4PhysicsVector::GetVectorLength() const
108 {
109  return numberOfNodes;
110 }
111 
112 //---------------------------------------------------------------
113 
114 inline
115  G4double G4PhysicsVector::GetValue(G4double theEnergy, G4bool&) const
116 {
117  size_t idx=0;
118  return Value(theEnergy, idx);
119 }
120 
121 //------------------------------------------------
122 
123 inline
124  G4double G4PhysicsVector::LinearInterpolation(size_t idx, G4double e) const
125 {
126  // Linear interpolation is used to get the value. Before this method
127  // is called it is ensured that the energy is inside the bin
128  // 0 < idx < numberOfNodes-1
129 
130  return dataVector[idx] +
131  ( dataVector[idx + 1]-dataVector[idx] ) * (e - binVector[idx])
132  /( binVector[idx + 1]-binVector[idx] );
133 }
134 
135 //---------------------------------------------------------------
136 
137 inline
138  G4double G4PhysicsVector::SplineInterpolation(size_t idx, G4double e) const
139 {
140  // Spline interpolation is used to get the value. Before this method
141  // is called it is ensured that the energy is inside the bin
142  // 0 < idx < numberOfNodes-1
143 
144  static const G4double onesixth = 1.0/6.0;
145 
146  // check bin value
147  G4double x1 = binVector[idx];
148  G4double x2 = binVector[idx + 1];
149  G4double delta = x2 - x1;
150 
151  G4double a = (x2 - e)/delta;
152  G4double b = (e - x1)/delta;
153 
154  // Final evaluation of cubic spline polynomial for return
155  G4double y1 = dataVector[idx];
156  G4double y2 = dataVector[idx + 1];
157 
158  G4double res = a*y1 + b*y2 +
159  ( (a*a*a - a)*secDerivative[idx] +
160  (b*b*b - b)*secDerivative[idx + 1] )*delta*delta*onesixth;
161 
162  return res;
163 }
164 
165 //---------------------------------------------------------------
166 
167 inline
168  G4double G4PhysicsVector::Interpolation(size_t idx, G4double e) const
169 {
170  G4double res;
171  if(useSpline) { res = SplineInterpolation(idx, e); }
172  else { res = LinearInterpolation(idx, e); }
173  return res;
174 }
175 
176 //---------------------------------------------------------------
177 
178 inline
179  void G4PhysicsVector::PutValue(size_t binNumber, G4double theValue)
180 {
181  dataVector[binNumber] = theValue;
182 }
183 
184 //---------------------------------------------------------------
185 
186 inline
187  G4bool G4PhysicsVector::IsFilledVectorExist() const
188 {
189  G4bool status=false;
190 
191  if(numberOfNodes > 0) { status=true; }
192  return status;
193 }
194 
195 //---------------------------------------------------------------
196 
197 inline
198  G4PhysicsVectorType G4PhysicsVector::GetType() const
199 {
200  return type;
201 }
202 
203 //---------------------------------------------------------------
204 
205 // Flag useSpline is "true" only if second derivatives are filled
206 inline
207  void G4PhysicsVector::SetSpline(G4bool val)
208 {
209  if(val) {
210  if(0 == secDerivative.size() && 0 < dataVector.size()) {
211  FillSecondDerivatives();
212  }
213  } else {
214  useSpline = false;
215  secDerivative.clear();
216  }
217 }
218 
219 //---------------------------------------------------------------
220 
221 inline
222 void G4PhysicsVector::SetVerboseLevel(G4int value)
223 {
224  verboseLevel = value;
225 }
226 
227 //---------------------------------------------------------------
228 
229 inline
230 G4int G4PhysicsVector::GetVerboseLevel(G4int)
231 {
232  return verboseLevel;
233 }
234 
235 //---------------------------------------------------------------
236 
237 inline
238 size_t G4PhysicsVector::FindBinLocation(G4double theEnergy) const
239 {
240  size_t bin;
241  if(type == T_G4PhysicsLogVector) {
242  bin = size_t(G4Log(theEnergy)/dBin - baseBin);
243  if(bin + 2 > numberOfNodes) { bin = numberOfNodes - 2; }
244  else if(bin > 0 && theEnergy < binVector[bin]) { --bin; }
245  else if(bin + 2 < numberOfNodes && theEnergy > binVector[bin+1])
246  { ++bin; }
247  } else if(type == T_G4PhysicsLinearVector) {
248  bin = size_t( theEnergy/dBin - baseBin );
249  if(bin + 2 > numberOfNodes) { bin = numberOfNodes - 2; }
250  else if(bin > 0 && theEnergy < binVector[bin]) { --bin; }
251  else if(bin + 2 < numberOfNodes && theEnergy > binVector[bin+1])
252  { ++bin; }
253  } else {
254  bin = 0;
255  size_t bin2;
256  size_t bin3 = numberOfNodes - 1;
257  while (bin != bin3 - 1) {
258  bin2 = bin + (bin3 - bin + 1)/2;
259  if (theEnergy > binVector[bin2]) { bin = bin2; }
260  else { bin3 = bin2; }
261  }
262 /*
263 // Bin location proposed by K.Genser (FNAL)
264 // V.I. Usage of this algorithm provides identical results
265 // no CPU advantage is observed in EM tests
266 // if new validation information will be known this code may be used
267  G4PVDataVector::const_iterator it =
268  std::lower_bound(binVector.begin(), binVector.end(), theEnergy);
269  bin = it - binVector.begin() - 1;
270 */
271  }
272  return bin;
273 }
274 
275 //---------------------------------------------------------------
276 
277 inline size_t G4PhysicsVector::FindBin(G4double e, size_t idx) const
278 {
279  size_t id = idx;
280  if(e < binVector[1]) {
281  id = 0;
282  } else if(e >= binVector[numberOfNodes-2]) {
283  id = numberOfNodes - 2;
284  } else if(idx >= numberOfNodes || e < binVector[idx] || e > binVector[idx+1]) {
285  id = FindBinLocation(e);
286  }
287  return id;
288 }
289 
290 //---------------------------------------------------------------