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

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