source: trunk/source/materials/src/G4IonStoppingData.cc@ 1353

Last change on this file since 1353 was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 15.9 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// $Id: G4IonStoppingData.cc,v 1.3 2010/10/25 08:41:39 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
28//
29// ===========================================================================
30// GEANT4 class source file
31//
32// Class: G4IonStoppingData
33//
34// Base class: G4VIonDEDXTable
35//
36// Author: Anton Lechner (Anton.Lechner@cern.ch)
37//
38// First implementation: 03. 11. 2009
39//
40// Modifications:
41// 25.10.2010 V.Ivanchenko fixed warnings reported by the Coverity tool
42//
43//
44// Class description: Class which can read ion stopping power data from
45// $G4LEDATA/ion_stopping_data
46//
47// Comments:
48//
49// ===========================================================================
50//
51
52#include "G4IonStoppingData.hh"
53#include "G4PhysicsVector.hh"
54#include "G4LPhysicsFreeVector.hh"
55#include <fstream>
56#include <sstream>
57#include <iomanip>
58
59
60// #########################################################################
61
62G4IonStoppingData::G4IonStoppingData(const G4String& leDirectory) :
63 subDir( leDirectory ) {
64
65}
66
67// #########################################################################
68
69G4IonStoppingData::~G4IonStoppingData() {
70
71 ClearTable();
72}
73
74// #########################################################################
75
76G4bool G4IonStoppingData::IsApplicable(
77 G4int atomicNumberIon, // Atomic number of ion
78 G4int atomicNumberElem // Atomic number of elemental material
79 ) {
80 G4bool isApplicable = true;
81 G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
82
83 G4IonDEDXMapElem::iterator iter = dedxMapElements.find(key);
84
85 if(iter == dedxMapElements.end()) isApplicable = false;
86
87 return isApplicable;
88}
89
90// #########################################################################
91
92G4bool G4IonStoppingData::IsApplicable(
93 G4int atomicNumberIon, // Atomic number of ion
94 const G4String& matIdentifier // Name or chemical formula of material
95 ) {
96 G4bool isApplicable = true;
97 G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
98
99 G4IonDEDXMapMat::iterator iter = dedxMapMaterials.find(key);
100
101 if(iter == dedxMapMaterials.end()) isApplicable = false;
102
103 return isApplicable;
104}
105
106// #########################################################################
107
108G4PhysicsVector* G4IonStoppingData::GetPhysicsVector(
109 G4int atomicNumberIon, // Atomic number of ion
110 G4int atomicNumberElem // Atomic number of elemental material
111 ) {
112
113 G4PhysicsVector* physVector = 0;
114
115 G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
116
117 G4IonDEDXMapElem::iterator iter = dedxMapElements.find(key);
118
119 if(iter != dedxMapElements.end()) physVector = iter -> second;
120
121 return physVector;
122}
123
124// #########################################################################
125
126G4PhysicsVector* G4IonStoppingData::GetPhysicsVector(
127 G4int atomicNumberIon, // Atomic number of ion
128 const G4String& matIdentifier // Name or chemical formula of material
129 ) {
130
131 G4PhysicsVector* physVector = 0;
132
133 G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
134
135 G4IonDEDXMapMat::iterator iter = dedxMapMaterials.find(key);
136
137 if(iter != dedxMapMaterials.end()) physVector = iter -> second;
138
139 return physVector;
140}
141
142// #########################################################################
143
144G4double G4IonStoppingData::GetDEDX(
145 G4double kinEnergyPerNucleon, // Kinetic energy per nucleon
146 G4int atomicNumberIon, // Atomic number of ion
147 G4int atomicNumberElem // Atomic number of elemental material
148 ) {
149 G4double dedx = 0;
150
151 G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
152
153 G4IonDEDXMapElem::iterator iter = dedxMapElements.find(key);
154
155 if( iter != dedxMapElements.end() ) {
156 G4PhysicsVector* physicsVector = iter -> second;
157
158 G4bool b;
159 dedx = physicsVector -> GetValue( kinEnergyPerNucleon, b );
160 }
161
162 return dedx;
163}
164
165// #########################################################################
166
167G4double G4IonStoppingData::GetDEDX(
168 G4double kinEnergyPerNucleon, // Kinetic energy per nucleon
169 G4int atomicNumberIon, // Atomic number of ion
170 const G4String& matIdentifier // Name or chemical formula of material
171 ) {
172 G4double dedx = 0;
173
174 G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
175
176 G4IonDEDXMapMat::iterator iter = dedxMapMaterials.find(key);
177
178 if(iter != dedxMapMaterials.end()) {
179 G4PhysicsVector* physicsVector = iter -> second;
180
181 G4bool b;
182 dedx = physicsVector -> GetValue( kinEnergyPerNucleon, b );
183 }
184
185 return dedx;
186}
187
188// #########################################################################
189
190G4bool G4IonStoppingData::AddPhysicsVector(
191 G4PhysicsVector* physicsVector, // Physics vector
192 G4int atomicNumberIon, // Atomic number of ion
193 const G4String& matIdentifier // Name of elemental material
194 ) {
195
196 if(physicsVector == 0) {
197
198#ifdef G4VERBOSE
199 G4cerr << "G4IonStoppingData::AddPhysicsVector() Error: Pointer to vector"
200 << " is null-pointer."
201 << G4endl;
202#endif
203
204 return false;
205 }
206
207 if(matIdentifier.empty()) {
208
209#ifdef G4VERBOSE
210 G4cerr << "G4IonStoppingData::AddPhysicsVector() Error: "
211 << "Cannot add physics vector. Invalid name."
212 << G4endl;
213#endif
214
215 return false;
216 }
217
218 if(atomicNumberIon <= 0) {
219
220#ifdef G4VERBOSE
221 G4cerr << "G4IonStoppingData::AddPhysicsVector() Error: "
222 << "Cannot add physics vector. Illegal atomic number."
223 << G4endl;
224#endif
225
226 return false;
227 }
228
229 G4IonDEDXKeyMat mkey = std::make_pair(atomicNumberIon, matIdentifier);
230
231 if(dedxMapMaterials.count(mkey) == 1) {
232
233#ifdef G4VERBOSE
234 G4cerr << "G4IonStoppingData::AddPhysicsVector() Error: "
235 << "Vector with Z1 = " << atomicNumberIon << ", mat = "
236 << matIdentifier
237 << "already exists. Remove first before replacing."
238 << G4endl;
239#endif
240
241 return false;
242 }
243
244 dedxMapMaterials[mkey] = physicsVector;
245
246 return true;
247}
248
249// #########################################################################
250
251G4bool G4IonStoppingData::AddPhysicsVector(
252 G4PhysicsVector* physicsVector, // Physics vector
253 G4int atomicNumberIon, // Atomic number of ion
254 G4int atomicNumberElem // Atomic number of elemental material
255 ) {
256
257 if(physicsVector == 0) {
258
259#ifdef G4VERBOSE
260 G4cerr << "G4IonStoppingData::AddPhysicsVector() Error: "
261 << "Pointer to vector is null-pointer."
262 << G4endl;
263#endif
264
265 return false;
266 }
267
268 if(atomicNumberIon <= 0) {
269
270#ifdef G4VERBOSE
271 G4cerr << "G4IonStoppingData::AddPhysicsVector() Error: "
272 << "Cannot add physics vector. Illegal atomic number."
273 << G4endl;
274#endif
275
276 return false;
277 }
278
279 if(atomicNumberElem <= 0) {
280
281#ifdef G4VERBOSE
282 G4cerr << "G4IonStoppingData::AddPhysicsVector() Error: "
283 << "Atomic number of element < 0."
284 << G4endl;
285#endif
286 return false;
287 }
288
289 G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
290
291 if(dedxMapElements.count(key) == 1) {
292
293#ifdef G4VERBOSE
294 G4cerr << "G4IonStoppingData::AddPhysicsVector() Error: "
295 << "Vector with Z1 = " << atomicNumberIon << ", Z2 = "
296 << atomicNumberElem
297 << " already exists. Remove first before replacing."
298 << G4endl;
299#endif
300 return false;
301 }
302
303 dedxMapElements[key] = physicsVector;
304
305 return true;
306}
307
308// #########################################################################
309
310G4bool G4IonStoppingData::RemovePhysicsVector(
311 G4int atomicNumberIon, // Atomic number of ion
312 const G4String& matIdentifier // Name or chemical formula of material
313 ) {
314
315 G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
316
317 G4IonDEDXMapMat::iterator iter = dedxMapMaterials.find(key);
318
319 if(iter == dedxMapMaterials.end()) {
320
321#ifdef G4VERBOSE
322 G4cerr << "G4IonStoppingData::RemovePhysicsVector() Warning: "
323 << "Cannot remove physics vector. Vector not found."
324 << G4endl;
325#endif
326
327 return false;
328 }
329
330 G4PhysicsVector* physicsVector = (*iter).second;
331
332 // Deleting key of physics vector from material map
333 dedxMapMaterials.erase(key);
334
335 // Deleting physics vector
336 delete physicsVector;
337
338 return true;
339}
340
341// #########################################################################
342
343G4bool G4IonStoppingData::RemovePhysicsVector(
344 G4int atomicNumberIon, // Atomic number of ion
345 G4int atomicNumberElem // Atomic number of elemental material
346 ) {
347 G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
348
349 G4IonDEDXMapElem::iterator iter = dedxMapElements.find(key);
350
351 if(iter == dedxMapElements.end()) {
352
353#ifdef G4VERBOSE
354 G4cerr << "G4IonStoppingData::RemovePhysicsVector() Warning: "
355 << "Cannot remove physics vector. Vector not found."
356 << G4endl;
357#endif
358
359 return false;
360 }
361
362 G4PhysicsVector* physicsVector = (*iter).second;
363
364 // Deleting key of physics vector from material map
365 dedxMapElements.erase(key);
366
367 // Deleting physics vector
368 delete physicsVector;
369
370 return true;
371}
372
373// #########################################################################
374
375G4bool G4IonStoppingData::BuildPhysicsVector(
376 G4int atomicNumberIon, // Atomic number of ion
377 const G4String& matIdentifier // Name of material
378 ) {
379
380 if( IsApplicable(atomicNumberIon, matIdentifier) ) return true;
381
382 char* path = getenv("G4LEDATA");
383 if ( !path ) {
384 G4Exception("G4IonStoppingData::BuildPhysicsVector: G4LEDATA environment variable not set");
385 return false;
386 }
387
388 std::ostringstream file;
389
390 file << path << "/" << subDir << "/z"
391 << atomicNumberIon << "_" << matIdentifier
392 << ".dat";
393
394 G4String fileName = G4String( file.str().c_str() );
395
396 std::ifstream ifilestream( fileName );
397
398 if ( !ifilestream.is_open() ) return false;
399
400 G4LPhysicsFreeVector* physicsVector = new G4LPhysicsFreeVector();
401
402 if( !physicsVector -> Retrieve(ifilestream, true) ) {
403
404 ifilestream.close();
405 return false;
406 }
407
408 physicsVector -> ScaleVector( MeV, MeV * cm2 /( 0.001 * g) );
409 physicsVector -> SetSpline( true );
410 physicsVector -> FillSecondDerivatives();
411
412 // Retrieved vector is added to material store
413 if( !AddPhysicsVector(physicsVector, atomicNumberIon, matIdentifier) ) {
414 delete physicsVector;
415 ifilestream.close();
416 return false;
417 }
418
419 ifilestream.close();
420 return true;
421}
422
423// #########################################################################
424
425G4bool G4IonStoppingData::BuildPhysicsVector(
426 G4int atomicNumberIon, // Atomic number of ion
427 G4int atomicNumberElem // Atomic number of elemental material
428 ) {
429
430 if( IsApplicable(atomicNumberIon, atomicNumberElem) ) return true;
431
432 char* path = getenv("G4LEDATA");
433 if ( !path ) {
434 G4Exception("G4IonStoppingData::BuildPhysicsVector: G4LEDATA environment variable not set");
435 return false;
436 }
437 std::ostringstream file;
438
439 file << path << "/" << subDir << "/z"
440 << atomicNumberIon << "_" << atomicNumberElem
441 << ".dat";
442
443 G4String fileName = G4String( file.str().c_str() );
444
445 std::ifstream ifilestream( fileName );
446
447 if ( !ifilestream.is_open() ) return false;
448
449 G4LPhysicsFreeVector* physicsVector = new G4LPhysicsFreeVector();
450
451 if( !physicsVector -> Retrieve(ifilestream, true) ) {
452
453 ifilestream.close();
454 return false;
455 }
456
457 physicsVector -> ScaleVector( MeV, MeV * cm2 /( 0.001 * g) );
458 physicsVector -> SetSpline( true );
459 physicsVector -> FillSecondDerivatives();
460
461 // Retrieved vector is added to material store
462 if( !AddPhysicsVector(physicsVector, atomicNumberIon, atomicNumberElem) ) {
463 delete physicsVector;
464 ifilestream.close();
465 return false;
466 }
467
468 ifilestream.close();
469 return true;
470}
471
472// #########################################################################
473
474void G4IonStoppingData::ClearTable() {
475
476 G4IonDEDXMapMat::iterator iterMat = dedxMapMaterials.begin();
477 G4IonDEDXMapMat::iterator iterMat_end = dedxMapMaterials.end();
478
479 for(;iterMat != iterMat_end; iterMat++) {
480
481 G4PhysicsVector* vec = iterMat -> second;
482
483 if(vec != 0) delete vec;
484 }
485
486 dedxMapMaterials.clear();
487
488 G4IonDEDXMapElem::iterator iterElem = dedxMapElements.begin();
489 G4IonDEDXMapElem::iterator iterElem_end = dedxMapElements.end();
490
491 for(;iterElem != iterElem_end; iterElem++) {
492
493 G4PhysicsVector* vec = iterElem -> second;
494
495 if(vec != 0) delete vec;
496 }
497
498 dedxMapElements.clear();
499}
500
501// #########################################################################
502
503void G4IonStoppingData::DumpMap() {
504
505 G4IonDEDXMapMat::iterator iterMat = dedxMapMaterials.begin();
506 G4IonDEDXMapMat::iterator iterMat_end = dedxMapMaterials.end();
507
508 G4cout << std::setw(15) << std::right
509 << "Atomic nmb ion"
510 << std::setw(25) << std::right
511 << "Material name"
512 << G4endl;
513
514 for(;iterMat != iterMat_end; iterMat++) {
515 G4IonDEDXKeyMat key = iterMat -> first;
516 G4PhysicsVector* physicsVector = iterMat -> second;
517
518 G4int atomicNumberIon = key.first;
519 G4String matIdentifier = key.second;
520
521 if(physicsVector != 0) {
522 G4cout << std::setw(15) << std::right
523 << atomicNumberIon
524 << std::setw(25) << std::right
525 << matIdentifier
526 << G4endl;
527 }
528 }
529
530 G4IonDEDXMapElem::iterator iterElem = dedxMapElements.begin();
531 G4IonDEDXMapElem::iterator iterElem_end = dedxMapElements.end();
532
533 G4cout << std::setw(15) << std::right
534 << "Atomic nmb ion"
535 << std::setw(25) << std::right
536 << "Atomic nmb material"
537 << G4endl;
538
539 for(;iterElem != iterElem_end; iterElem++) {
540 G4IonDEDXKeyElem key = iterElem -> first;
541 G4PhysicsVector* physicsVector = iterElem -> second;
542
543 G4int atomicNumberIon = key.first;
544 G4int atomicNumberElem = key.second;
545
546 if(physicsVector != 0) {
547 G4cout << std::setw(15) << std::right
548 << atomicNumberIon
549 << std::setw(25) << std::right
550 << atomicNumberElem
551 << G4endl;
552 }
553 }
554
555}
556
557// #########################################################################
558
Note: See TracBrowser for help on using the repository browser.