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

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