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

Last change on this file since 1182 was 1058, checked in by garnier, 15 years ago

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