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

Last change on this file since 1302 was 1196, checked in by garnier, 16 years ago

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