source: trunk/source/materials/src/G4MaterialPropertyVector.cc@ 896

Last change on this file since 896 was 850, checked in by garnier, 17 years ago

geant4.8.2 beta
File size: 11.3 KB
Line 
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: G4MaterialPropertyVector.cc,v 1.16 2008/06/05 23:38:51 gum Exp $
28// GEANT4 tag $Name: HEAD $
29//
30//
31////////////////////////////////////////////////////////////////////////
32// G4MaterialPropertyVector Class Implementation
33////////////////////////////////////////////////////////////////////////
34//
35// File: G4MaterialPropertyVector.cc
36// Version: 1.0
37// Created: 1996-02-08
38// Author: Juliet Armstrong
39// Updated: 1997-03-25 by Peter Gumplinger
40// > cosmetics (only)
41// mail: gum@triumf.ca
42//
43////////////////////////////////////////////////////////////////////////
44
45#include "G4MaterialPropertyVector.hh"
46
47/////////////////////////
48// Class Implementation
49/////////////////////////
50
51 //////////////
52 // Operators
53 //////////////
54
55// ++ operator
56// -----------
57//
58G4bool G4MaterialPropertyVector::operator ++()
59{
60 CurrentEntry++;
61 if (CurrentEntry < NumEntries) {
62 return true;
63 }
64 else {
65 return false;
66 }
67}
68
69// = operator
70// ----------
71//
72G4MaterialPropertyVector&
73G4MaterialPropertyVector::operator =(const G4MaterialPropertyVector& right)
74{
75 if (this == &right) return *this;
76
77 // clear the vector of current contents
78
79 MPV.clear();
80
81 // create an actual copy (instead of the shallow copy that the
82 // assignment operator defaults to for G4RWTPtrSortedVector)
83
84 NumEntries = 0;
85 CurrentEntry = -1;
86
87 for (G4int i = 0 ; i < right.NumEntries; i++) {
88 G4MPVEntry *newElement = new G4MPVEntry(right.GetEntry(i));
89 MPV.push_back(newElement);
90 NumEntries++;
91 }
92
93 return *this;
94}
95
96 /////////////////
97 // Constructors
98 /////////////////
99
100G4MaterialPropertyVector::G4MaterialPropertyVector(G4double *PhotonEnergies,
101 G4double *PropertyValues,
102 G4int NumElements)
103{
104 NumEntries = 0;
105 CurrentEntry = -1;
106
107 // create a vector filling it with the values
108 // from PhotonEnergies[] and PropertyValues[]
109
110 for(G4int i = 0; i < NumElements; i++) {
111 AddElement(PhotonEnergies[i], PropertyValues[i]);
112 }
113}
114
115G4MaterialPropertyVector::G4MaterialPropertyVector
116 (const G4MaterialPropertyVector &right)
117{
118 // create an actual copy (instead of the shallow copy that the
119 // assignment operator defaults to for G4RWTPtrSortedVector)
120
121 NumEntries = 0;
122 CurrentEntry = -1;
123
124 for (G4int i = 0 ; i < right.NumEntries; i++) {
125 G4MPVEntry *newElement = new G4MPVEntry(right.GetEntry(i));
126 MPV.push_back(newElement);
127 NumEntries++;
128 }
129}
130
131 ////////////////
132 // Destructors
133 ////////////////
134
135G4MaterialPropertyVector::~G4MaterialPropertyVector()
136{
137 MPV.clear();
138}
139
140 ////////////
141 // Methods
142 ////////////
143void G4MaterialPropertyVector::AddElement(G4double aPhotonEnergy,
144 G4double aPropertyValue)
145{
146 G4MPVEntry *newElement;
147
148 newElement = new G4MPVEntry(aPhotonEnergy, aPropertyValue);
149 MPV.push_back(newElement);
150 std::sort(MPV.begin(), MPV.end(), MPVEntry_less());
151 NumEntries++;
152}
153
154void G4MaterialPropertyVector::RemoveElement(G4double aPhotonEnergy)
155{
156 G4MPVEntry *newElement;
157 G4MPVEntry *success=0;
158
159 newElement = new G4MPVEntry(aPhotonEnergy, DBL_MAX);
160
161 std::vector<G4MPVEntry*>::iterator i;
162 for (i = MPV.begin(); i != MPV.end(); i++)
163 if (**i == *newElement) {success = *i; break;}
164 // success = MPV.remove(newElement);
165
166 if(success == 0)
167 {
168 G4Exception("G4MaterialPropertyVector::RemoveElement==>"
169 "element not found");
170 return;
171 }
172 else
173 {
174 MPV.erase(i); // remove done here.
175 }
176
177 NumEntries--;
178}
179
180G4double
181G4MaterialPropertyVector::GetProperty(G4double aPhotonEnergy) const
182{
183 G4MPVEntry *target, *temp;
184 G4int left, right;
185 G4double ratio1, ratio2, pmright, pmleft, InterpolatedValue;
186
187 /////////////////////////
188 // Establish table range
189 /////////////////////////
190
191 G4double PMmin = MPV.front()->GetPhotonEnergy();
192 G4double minProp = MPV.front()->GetProperty();
193 G4double PMmax = MPV.back() ->GetPhotonEnergy();
194 G4double maxProp = MPV.back() ->GetProperty();
195
196 ///////////////////////////////////////////
197 // Does value fall outside range of table?
198 ///////////////////////////////////////////
199
200 if (aPhotonEnergy < PMmin)
201 {
202 G4cout << "\nWarning: G4MaterialPropertyVector::GetProperty";
203 G4cout << "\n==> attempt to Retrieve Property below range"
204 << G4endl;
205 return minProp;
206 }
207
208 if (aPhotonEnergy > PMmax)
209 {
210 G4cout << "\nWarning: G4MaterialPropertyVector::GetProperty";
211 G4cout << "\n==> attempt to Retrieve Property above range"
212 << G4endl;
213 return maxProp;
214 }
215
216 target = new G4MPVEntry(aPhotonEnergy, 0.0);
217
218 temp = 0;
219 //temp = MPV.find(target);
220 std::vector<G4MPVEntry*>::const_iterator i;
221 for (i = MPV.begin(); i != MPV.end(); i++)
222 if (**i == *target) {temp = *i; break;}
223 if (temp != 0) {
224
225 ////////////////////////
226 // Return actual value
227 ////////////////////////
228
229 G4double retval = temp->GetProperty();
230 delete target;
231 return retval;
232 }
233 else {
234 //////////////////////////////
235 // Return interpolated value
236 //////////////////////////////
237
238 GetAdjacentBins(aPhotonEnergy, &left, &right);
239
240 pmleft = MPV[left]->GetPhotonEnergy();
241 pmright = MPV[right]->GetPhotonEnergy();
242 ratio1 = (aPhotonEnergy-pmleft)/(pmright-pmleft);
243 ratio2 = 1 - ratio1;
244 InterpolatedValue = MPV[left]->GetProperty()*ratio2 +
245 MPV[right]->GetProperty()*ratio1;
246
247 delete target;
248 return InterpolatedValue;
249 }
250}
251
252G4double G4MaterialPropertyVector::GetProperty() const
253{
254// For use with G4MaterialPropertyVector iterator
255
256 if(CurrentEntry == -1 || CurrentEntry >= NumEntries) {
257 G4Exception("G4MaterialPropertyVector::GetProperty ==>"
258 "Iterator attempted to Retrieve Property out of range");
259 return DBL_MAX;
260 }
261 else {
262 return MPV[CurrentEntry]->GetProperty();
263 }
264}
265
266G4double G4MaterialPropertyVector::GetPhotonEnergy() const
267{
268// For use with G4MaterialPropertyVector iterator
269
270 if(CurrentEntry == -1 || CurrentEntry >= NumEntries) {
271 G4Exception("G4MaterialPropertyVector::GetPhotonEnergy ==>"
272 "Iterator attempted to Retrieve Photon Energy out of range");
273 return DBL_MAX;
274 }
275 else {
276 return MPV[CurrentEntry]->GetPhotonEnergy();
277 }
278}
279
280G4double
281G4MaterialPropertyVector::GetPhotonEnergy(G4double aProperty) const
282{
283// ***NB***
284// Assumes that the property is an increasing function of photon energy (e.g.
285// refraction index)
286// ***NB***
287//
288// Returns the photon energy corresponding to the property value passed in.
289// If several photon energy values correspond to the value passed in, the
290// function returns the first photon energy in the vector that corresponds
291// to that value.
292
293 G4int left, right, mid;
294 G4double ratio1, ratio2, pright, pleft, InterpolatedValue;
295
296 //////////////////////////
297 // Establish Table range
298 //////////////////////////
299
300 G4double PropMin = MPV.front()->GetProperty();
301 G4double PMmin = MPV.front()->GetPhotonEnergy();
302 G4double PropMax = MPV.back() ->GetProperty();
303 G4double PMmax = MPV.back() ->GetPhotonEnergy();
304
305 ///////////////////////////////////////////
306 // Does value fall outside range of table?
307 ///////////////////////////////////////////
308
309 if (aProperty < PropMin)
310 {
311 G4cout << "\nWarning: G4MaterialPropertyVector::GetPhotonEnergy";
312 G4cout << "\n==> attempt to Retrieve Photon Energy out of range"
313 << G4endl;
314 return PMmin;
315 }
316
317 if (aProperty > PropMax) {
318 G4cout << "\nWarning: G4MaterialPropertyVector::GetPhotonEnergy";
319 G4cout << "\n==> attempt to Retrieve Photon Energy out of range"
320 << G4endl;
321 return PMmax;
322 }
323
324 //////////////////////////////
325 // Return interpolated value
326 //////////////////////////////
327
328 left = 0;
329 right = MPV.size(); // was .entries()
330
331 // find values in bins on either side of aProperty
332
333 do {
334 mid = (left + right)/2;
335 if (MPV[mid]->GetProperty() == aProperty) {
336
337 // Get first photon energy value in vector that
338 // corresponds to property value
339
340 while (mid-1 >= 0 &&
341 MPV[mid-1]->GetProperty() == aProperty) {
342 mid--;
343 }
344
345 InterpolatedValue = MPV[mid]->GetPhotonEnergy();
346 goto end_GetPhotonEnergy;
347 }
348 if (MPV[mid]->GetProperty() < aProperty)
349 left = mid;
350 else
351 right = mid;
352
353 } while ((right - left) > 1);
354
355 pleft = MPV[left]->GetProperty();
356 pright = MPV[right]->GetProperty();
357 ratio1 = (aProperty - pleft) / (pright - pleft);
358 ratio2 = 1 - ratio1;
359 InterpolatedValue = MPV[left]->GetPhotonEnergy()*ratio2 +
360 MPV[right]->GetPhotonEnergy()*ratio1;
361
362end_GetPhotonEnergy:
363
364 return InterpolatedValue;
365}
366
367void G4MaterialPropertyVector::ResetIterator()
368{
369 CurrentEntry = -1;
370}
371
372void G4MaterialPropertyVector::DumpVector()
373{
374 if (MPV.empty()) {
375 G4cerr << "nothing to dump" << G4endl;
376 G4Exception("G4MaterialPropertyVector::DumpVector ==>"
377 "Nothing to dump! Vector is empty");
378 }
379
380 for (G4int i = 0; i < NumEntries; i++) {
381 G4cout << "MPV["<< i << "]: ";
382 MPV[i]->DumpEntry();
383 }
384 G4cout << " Done DumpVector of " << NumEntries << " entries " << G4endl;
385
386}
387
388G4MPVEntry G4MaterialPropertyVector::GetEntry(G4int i) const
389{
390 return *MPV[i];
391}
392
393void
394G4MaterialPropertyVector::GetAdjacentBins(G4double aPhotonEnergy,
395 G4int *left,
396 G4int *right) const
397{
398 G4int mid;
399
400 *left = 0;
401 *right = (MPV.size() - 1); // was .entries()
402
403 // find values in bins on either side of aPhotonEnergy
404
405 do {
406 mid = (*left + *right)/2;
407
408 if (MPV[mid]->GetPhotonEnergy() < aPhotonEnergy)
409 {
410 *left = mid;
411 }
412 else
413 {
414 *right = mid;
415 }
416 } while ((*right - *left) > 1);
417}
Note: See TracBrowser for help on using the repository browser.