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

Last change on this file since 1347 was 1347, checked in by garnier, 13 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.