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

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

geant4 tag 9.4
File size: 17.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: G4ExtDEDXTable.cc,v 1.4 2010/11/01 18:18:57 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
28//
29// ===========================================================================
30// GEANT4 class source file
31//
32// Class:                G4ExtDEDXTable
33//
34// Base class:           G4VIonDEDXTable
35//
36// Author:               Anton Lechner (Anton.Lechner@cern.ch)
37//
38// First implementation: 29. 02. 2009
39//
40// Modifications:
41// 03.11.2009 A. Lechner:  Added new methods BuildPhysicsVector according
42//            to interface changes in base class G4VIonDEDXTable.
43// 25.10.2010 V.Ivanchenko fixed bug in usage of iterators reported by the
44//            Coverity tool
45// 01.11.2010 V.Ivanchenko fixed remaining bugs reported by Coverity
46//
47//
48// Class description:
49//    Utility class for users to add their own electronic stopping powers
50//    for ions. This class is dedicated for use with G4IonParametrisedLossModel
51//    of the low-energy electromagnetic package.
52//
53// Comments:
54//
55// ===========================================================================
56//
57
58#include "G4ExtDEDXTable.hh"
59#include "G4PhysicsVector.hh"
60#include "G4PhysicsVectorType.hh"
61#include "G4LPhysicsFreeVector.hh"
62#include "G4PhysicsLogVector.hh"
63#include "G4PhysicsFreeVector.hh"
64#include "G4PhysicsOrderedFreeVector.hh"
65#include "G4PhysicsLinearVector.hh"
66#include "G4PhysicsLnVector.hh"
67#include <fstream>
68#include <sstream>
69#include <iomanip>
70
71
72// #########################################################################
73
74G4ExtDEDXTable::G4ExtDEDXTable() {
75
76}
77
78// #########################################################################
79
80G4ExtDEDXTable::~G4ExtDEDXTable() {
81
82  ClearTable();
83}
84
85// #########################################################################
86
87G4bool G4ExtDEDXTable::BuildPhysicsVector(G4int ionZ, G4int matZ) {
88
89  return IsApplicable( ionZ, matZ );
90}
91
92
93// #########################################################################
94
95G4bool G4ExtDEDXTable::BuildPhysicsVector(G4int ionZ, 
96                                          const G4String& matName) {
97
98  return IsApplicable( ionZ, matName );
99}
100
101// #########################################################################
102
103G4bool G4ExtDEDXTable::IsApplicable(
104         G4int atomicNumberIon,  // Atomic number of ion
105         G4int atomicNumberElem  // Atomic number of elemental material
106                                    ) {
107  G4bool isApplicable = true; 
108  G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
109
110  G4IonDEDXMapElem::iterator iter = dedxMapElements.find(key);
111
112  if(iter == dedxMapElements.end()) isApplicable = false; 
113
114  return isApplicable; 
115}
116
117// #########################################################################
118
119G4bool G4ExtDEDXTable::IsApplicable(
120         G4int atomicNumberIon,         // Atomic number of ion
121         const G4String& matIdentifier  // Name or chemical formula of material
122                                    ) {
123  G4bool isApplicable = true; 
124  G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
125
126  G4IonDEDXMapMat::iterator iter = dedxMapMaterials.find(key);
127
128  if(iter == dedxMapMaterials.end()) isApplicable = false; 
129
130  return isApplicable; 
131}
132
133// #########################################################################
134
135G4PhysicsVector* G4ExtDEDXTable::GetPhysicsVector(
136         G4int atomicNumberIon,        // Atomic number of ion
137         G4int atomicNumberElem        // Atomic number of elemental material
138                                    ) {
139
140  G4PhysicsVector* physVector = 0;
141
142  G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
143
144  G4IonDEDXMapElem::iterator iter = dedxMapElements.find(key);
145
146  if(iter != dedxMapElements.end()) physVector = iter -> second; 
147
148  return physVector; 
149}
150
151// #########################################################################
152
153G4PhysicsVector*  G4ExtDEDXTable::GetPhysicsVector(
154         G4int atomicNumberIon,        // Atomic number of ion
155         const G4String& matIdentifier // Name or chemical formula of material
156                                    ) {
157
158  G4PhysicsVector* physVector = 0;
159
160  G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
161
162  G4IonDEDXMapMat::iterator iter = dedxMapMaterials.find(key);
163
164  if(iter != dedxMapMaterials.end()) physVector = iter -> second; 
165
166  return physVector; 
167}
168
169// #########################################################################
170
171G4double G4ExtDEDXTable::GetDEDX(
172         G4double kinEnergyPerNucleon, // Kinetic energy per nucleon
173         G4int atomicNumberIon,        // Atomic number of ion
174         G4int atomicNumberElem        // Atomic number of elemental material
175                                  ) {
176  G4double dedx = 0;
177
178  G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
179
180  G4IonDEDXMapElem::iterator iter = dedxMapElements.find(key);
181
182  if( iter != dedxMapElements.end() ) {
183     G4PhysicsVector* physicsVector = iter -> second; 
184
185     G4bool b;
186     dedx = physicsVector -> GetValue( kinEnergyPerNucleon, b );   
187  }
188
189  return dedx; 
190}
191
192// #########################################################################
193
194G4double G4ExtDEDXTable::GetDEDX(
195         G4double kinEnergyPerNucleon, // Kinetic energy per nucleon
196         G4int atomicNumberIon,        // Atomic number of ion
197         const G4String& matIdentifier // Name or chemical formula of material
198                                  ) {
199  G4double dedx = 0;
200
201  G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
202
203  G4IonDEDXMapMat::iterator iter = dedxMapMaterials.find(key);
204
205  if(iter != dedxMapMaterials.end()) {
206     G4PhysicsVector* physicsVector = iter -> second; 
207
208     G4bool b;
209     dedx = physicsVector -> GetValue( kinEnergyPerNucleon, b );   
210  }
211
212  return dedx; 
213}
214
215// #########################################################################
216
217G4bool G4ExtDEDXTable::AddPhysicsVector(
218        G4PhysicsVector* physicsVector, // Physics vector
219        G4int atomicNumberIon,          // Atomic number of ion
220        const G4String& matIdentifier,  // Name of elemental material
221        G4int atomicNumberElem          // Atomic number of elemental material
222                                      ) {
223
224  if(physicsVector == 0) {
225
226#ifdef G4VERBOSE
227     G4cout << "G4IonDEDXTable::AddPhysicsVector() Error: Pointer to vector"
228            << " is null-pointer."
229            << G4endl;
230#endif
231
232     return false;
233  }
234
235  if(matIdentifier.empty()) {
236
237#ifdef G4VERBOSE
238     G4cout << "G4IonDEDXTable::AddPhysicsVector() Error: "
239            << "Cannot add physics vector. Invalid name."
240            << G4endl;
241#endif
242
243     return false;
244  }
245
246  if(atomicNumberIon <= 2) {
247
248#ifdef G4VERBOSE
249     G4cout << "G4IonDEDXTable::AddPhysicsVector() Error: "
250            << "Cannot add physics vector. Illegal atomic number."
251            << G4endl;
252#endif
253
254     return false;
255  }
256
257  if(atomicNumberElem > 0) {
258
259     G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
260
261     if(dedxMapElements.count(key) == 1) {
262
263#ifdef G4VERBOSE
264        G4cout << "G4IonDEDXTable::AddPhysicsVector() Error: "
265               << "Vector already exists. Remove first before replacing."
266               << G4endl;
267#endif
268        return false;
269     }
270
271     dedxMapElements[key] = physicsVector;
272  }
273
274  G4IonDEDXKeyMat mkey = std::make_pair(atomicNumberIon, matIdentifier);
275
276  if(dedxMapMaterials.count(mkey) == 1) {
277
278#ifdef G4VERBOSE
279     G4cout << "G4IonDEDXTable::AddPhysicsVector() Error: "
280            << "Vector already exists. Remove first before replacing."
281            << G4endl;
282#endif
283
284     return false;
285  }
286
287  dedxMapMaterials[mkey] = physicsVector;
288
289  return true;
290}
291
292// #########################################################################
293
294G4bool G4ExtDEDXTable::RemovePhysicsVector(
295        G4int atomicNumberIon,         // Atomic number of ion
296        const G4String& matIdentifier  // Name or chemical formula of material
297                                      ) {
298
299  G4PhysicsVector* physicsVector = 0;
300
301  // Deleting key of physics vector from material map
302  G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
303
304  G4IonDEDXMapMat::iterator iter = dedxMapMaterials.find(key);
305
306  if(iter == dedxMapMaterials.end()) {
307
308#ifdef G4VERBOSE
309    G4cout << "G4IonDEDXTable::RemovePhysicsVector() Warning: "
310           << "Cannot remove physics vector. Vector not found."
311           << G4endl;
312#endif
313
314    return false;
315  }
316
317  physicsVector = (*iter).second;
318  dedxMapMaterials.erase(key);
319
320  // Deleting key of physics vector from elemental material map (if it exists)
321  G4IonDEDXMapElem::iterator it;
322 
323  for(it=dedxMapElements.begin(); it!=dedxMapElements.end(); ++it) {
324
325     if( (*it).second == physicsVector ) {
326        dedxMapElements.erase(it);
327        break;
328     }
329  }
330
331  // Deleting physics vector
332  delete physicsVector;
333
334  return true;
335}
336
337// #########################################################################
338
339G4bool G4ExtDEDXTable::StorePhysicsTable(
340         const G4String& fileName // File name
341                                    ) {
342  G4bool success = true;
343
344  std::ofstream ofilestream;
345
346  ofilestream.open( fileName, std::ios::out );
347
348  if( !ofilestream ) {
349
350#ifdef G4VERBOSE
351     G4cout << "G4ExtDEDXTable::StorePhysicsVector() " 
352            << " Cannot open file "<< fileName
353            << G4endl;
354#endif
355     
356     success = false;
357  }   
358  else {
359
360     size_t nmbMatTables = dedxMapMaterials.size();
361
362     ofilestream << nmbMatTables << G4endl << G4endl; 
363
364     G4IonDEDXMapMat::iterator iterMat = dedxMapMaterials.begin();
365     G4IonDEDXMapMat::iterator iterMat_end = dedxMapMaterials.end();
366
367     for(;iterMat != iterMat_end; iterMat++) {
368         G4IonDEDXKeyMat key = iterMat -> first;
369         G4PhysicsVector* physicsVector = iterMat -> second; 
370
371         G4int atomicNumberIon = key.first;
372         G4String matIdentifier = key.second;
373
374         G4int atomicNumberElem = FindAtomicNumberElement(physicsVector);
375
376         if(physicsVector != 0) {
377            ofilestream << atomicNumberIon << "  " << matIdentifier;
378
379            if(atomicNumberElem > 0) ofilestream << "  " << atomicNumberElem;
380
381            ofilestream << "  # <Atomic number ion>  <Material name>  ";
382
383            if(atomicNumberElem > 0) ofilestream << "<Atomic number element>";
384
385            ofilestream << G4endl << physicsVector -> GetType() << G4endl;
386
387            physicsVector -> Store(ofilestream, true);
388
389            ofilestream << G4endl;
390         }
391         else {
392
393#ifdef G4VERBOSE
394              G4cout << "G4ExtDEDXTable::StorePhysicsVector() " 
395                     << " Cannot store physics vector." 
396                     << G4endl;
397#endif
398
399         }
400     }
401  }
402
403  ofilestream.close();
404
405  return success; 
406}
407
408// #########################################################################
409
410G4bool G4ExtDEDXTable::RetrievePhysicsTable(
411              const G4String& fileName  // File name
412                                             ) { 
413
414  std::ifstream ifilestream;
415  ifilestream.open( fileName, std::ios::in );
416  if( ! ifilestream ) {
417#ifdef G4VERBOSE
418    G4cout << "G4ExtDEDXTable::RetrievePhysicsVector() " 
419           << " Cannot open file "<< fileName
420           << G4endl;
421#endif
422    return false;
423  }   
424
425  std::string::size_type nmbVectors = 0;
426  ifilestream >> nmbVectors;
427
428  if(nmbVectors == std::string::npos || nmbVectors == 0) {
429#ifdef G4VERBOSE
430    G4cout << "G4ExtDEDXTable::RetrievePhysicsVector() " 
431           << " The file is empty "<< nmbVectors
432           << G4endl;
433#endif
434    return false;
435  } 
436  for(size_t i = 0; i<nmbVectors; ++i) {
437
438    G4String line = "";
439    while( line.empty() ) {
440
441      getline( ifilestream, line );
442      if( ifilestream.fail() ) { 
443#ifdef G4VERBOSE 
444        G4cout << "G4ExtDEDXTable::RetrievePhysicsTable() " 
445               << " File content of " << fileName << " ill-formated." 
446               << G4endl;
447#endif         
448        ifilestream.close(); 
449        return false; 
450      }
451
452      std::string::size_type pos = line.find_first_of("#");
453      if(pos != std::string::npos && pos > 0) {
454        line = line.substr(0, pos);
455      }
456    }
457
458    std::istringstream headerstream( line );     
459
460    std::string::size_type atomicNumberIon;
461    headerstream >> atomicNumberIon;
462
463    G4String materialName;
464    headerstream >> materialName;
465
466    if( headerstream.fail() || std::string::npos == atomicNumberIon) {
467 
468#ifdef G4VERBOSE 
469      G4cout << "G4ExtDEDXTable::RetrievePhysicsTable() " 
470             << " File content of " << fileName << " ill-formated "
471             << " (vector header)." 
472             << G4endl;
473#endif         
474      ifilestream.close();
475      return false;
476    } 
477
478    std::string::size_type atomicNumberMat;
479    headerstream >> atomicNumberMat;
480
481    if( headerstream.eof() || std::string::npos == atomicNumberMat) { 
482      atomicNumberMat = 0; 
483    }
484
485    G4int vectorType;
486    ifilestream >> vectorType;
487     
488    G4PhysicsVector* physicsVector = CreatePhysicsVector(vectorType);
489
490    if(physicsVector == 0) {
491#ifdef G4VERBOSE 
492      G4cout << "G4ExtDEDXTable::RetrievePhysicsTable  "
493             << " illegal physics Vector type " << vectorType
494             << " in  " << fileName
495             << G4endl;
496#endif         
497      ifilestream.close();
498      return false;
499    }
500
501    if( !physicsVector -> Retrieve(ifilestream, true) ) {
502       
503#ifdef G4VERBOSE 
504      G4cout << "G4ExtDEDXTable::RetrievePhysicsTable() " 
505             << " File content of " << fileName << " ill-formated." 
506             << G4endl;
507#endif         
508      ifilestream.close();
509      return false;
510    } 
511
512    physicsVector -> SetSpline(true);
513
514    // Retrieved vector is added to material store
515    if( !AddPhysicsVector(physicsVector, (G4int)atomicNumberIon, 
516                          materialName, (G4int)atomicNumberMat) ) {
517
518      delete physicsVector;
519      ifilestream.close();
520      return false;
521    }
522  }
523
524  ifilestream.close();
525
526  return true;
527}
528
529// #########################################################################
530
531G4PhysicsVector* G4ExtDEDXTable::CreatePhysicsVector(G4int vectorType) {
532
533  G4PhysicsVector* physicsVector = 0;
534
535  switch (vectorType) {
536
537  case T_G4PhysicsLinearVector: 
538    physicsVector = new G4PhysicsLinearVector();
539    break;
540
541  case T_G4PhysicsLogVector: 
542    physicsVector = new G4PhysicsLogVector();
543    break;
544
545  case T_G4PhysicsLnVector: 
546    physicsVector = new G4PhysicsLnVector();
547    break;
548
549  case T_G4PhysicsFreeVector: 
550    physicsVector = new G4PhysicsFreeVector();
551    break;
552
553  case T_G4PhysicsOrderedFreeVector: 
554    physicsVector = new G4PhysicsOrderedFreeVector();
555    break;
556
557  case T_G4LPhysicsFreeVector: 
558    physicsVector = new G4LPhysicsFreeVector();
559    break;
560 
561  default:
562    break;
563  }
564  return physicsVector;
565}
566
567// #########################################################################
568
569G4int G4ExtDEDXTable::FindAtomicNumberElement(
570                                   G4PhysicsVector* physicsVector
571                                                   ) {
572
573  G4int atomicNumber = 0;
574
575  G4IonDEDXMapElem::iterator iter = dedxMapElements.begin();
576  G4IonDEDXMapElem::iterator iter_end = dedxMapElements.end();
577 
578  for(;iter != iter_end; iter++) {
579
580     if( (*iter).second == physicsVector ) {
581
582        G4IonDEDXKeyElem key = (*iter).first;
583        atomicNumber = key.second;
584     }
585  }
586
587  return atomicNumber;
588}
589
590// #########################################################################
591
592void G4ExtDEDXTable::ClearTable() {
593
594  G4IonDEDXMapMat::iterator iterMat = dedxMapMaterials.begin();
595  G4IonDEDXMapMat::iterator iterMat_end = dedxMapMaterials.end();
596
597  for(;iterMat != iterMat_end; iterMat++) { 
598
599    G4PhysicsVector* vec = iterMat -> second;
600
601    if(vec != 0) delete vec;
602  }
603
604  dedxMapElements.clear();
605  dedxMapMaterials.clear();
606}
607
608// #########################################################################
609
610void G4ExtDEDXTable::DumpMap() {
611
612  G4IonDEDXMapMat::iterator iterMat = dedxMapMaterials.begin();
613  G4IonDEDXMapMat::iterator iterMat_end = dedxMapMaterials.end();
614
615  G4cout << std::setw(15) << std::right
616         << "Atomic nmb ion"
617         << std::setw(25) << std::right
618         << "Material name"
619         << std::setw(25) << std::right
620         << "Atomic nmb material"
621         << G4endl;
622
623  for(;iterMat != iterMat_end; iterMat++) {
624      G4IonDEDXKeyMat key = iterMat -> first;
625      G4PhysicsVector* physicsVector = iterMat -> second; 
626
627      G4int atomicNumberIon = key.first;
628      G4String matIdentifier = key.second;
629
630      G4int atomicNumberElem = FindAtomicNumberElement(physicsVector);
631
632      if(physicsVector != 0) {
633         G4cout << std::setw(15) << std::right
634                << atomicNumberIon
635                << std::setw(25) << std::right
636                << matIdentifier
637                << std::setw(25) << std::right;
638
639         if(atomicNumberElem > 0) G4cout << atomicNumberElem;
640         else G4cout << "N/A";
641
642         G4cout << G4endl;
643      }
644  }
645
646}
647
648// #########################################################################
649
Note: See TracBrowser for help on using the repository browser.