| Geant4 Cross Reference |
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 // GEANT4 class source file
28 //
29 // Class: G4ExtDEDXTable
30 //
31 // Base class: G4VIonDEDXTable
32 //
33 // Author: Anton Lechner (Anton.Lechner@cern.ch)
34 //
35 // First implementation: 29. 02. 2009
36 //
37 // Modifications:
38 // 03.11.2009 A. Lechner: Added new methods BuildPhysicsVector according
39 // to interface changes in base class G4VIonDEDXTable.
40 // 25.10.2010 V.Ivanchenko fixed bug in usage of iterators reported by the
41 // Coverity tool
42 // 01.11.2010 V.Ivanchenko fixed remaining bugs reported by Coverity
43 //
44 //
45 // Class description:
46 // Utility class for users to add their own electronic stopping powers
47 // for ions. This class is dedicated for use with G4IonParametrisedLossModel
48 // of the low-energy electromagnetic package.
49 //
50 // Comments:
51 //
52 // ===========================================================================
53
54 #include "G4ExtDEDXTable.hh"
55
56 #include "G4PhysicsFreeVector.hh"
57 #include "G4PhysicsLinearVector.hh"
58 #include "G4PhysicsLogVector.hh"
59 #include "G4PhysicsVector.hh"
60 #include "G4PhysicsVectorType.hh"
61
62 #include <fstream>
63 #include <iomanip>
64 #include <sstream>
65
66 // #########################################################################
67
68 G4ExtDEDXTable::~G4ExtDEDXTable() { ClearTable(); }
69
70 // #########################################################################
71
72 G4bool G4ExtDEDXTable::BuildPhysicsVector(G4int ionZ, G4int matZ)
73 {
74 return IsApplicable(ionZ, matZ);
75 }
76
77 // #########################################################################
78
79 G4bool G4ExtDEDXTable::BuildPhysicsVector(G4int ionZ, const G4String& matName)
80 {
81 return IsApplicable(ionZ, matName);
82 }
83
84 // #########################################################################
85
86 G4bool G4ExtDEDXTable::IsApplicable(G4int atomicNumberIon, // Atomic number of ion
87 G4int atomicNumberElem // Atomic number of elemental material
88 )
89 {
90 G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
91
92 auto iter = dedxMapElements.find(key);
93
94 return iter != dedxMapElements.end();
95 }
96
97 // #########################################################################
98
99 G4bool G4ExtDEDXTable::IsApplicable(G4int atomicNumberIon, // Atomic number of ion
100 const G4String& matIdentifier // Name or chemical formula of material
101 )
102 {
103 G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
104
105 auto iter = dedxMapMaterials.find(key);
106
107 return iter != dedxMapMaterials.end();
108 }
109
110 // #########################################################################
111
112 G4PhysicsVector* G4ExtDEDXTable::GetPhysicsVector(G4int atomicNumberIon, // Atomic number of ion
113 G4int atomicNumberElem // Atomic number of elemental material
114 )
115 {
116 G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
117
118 auto iter = dedxMapElements.find(key);
119
120 return (iter != dedxMapElements.end()) ? iter->second : nullptr;
121 }
122
123 // #########################################################################
124
125 G4PhysicsVector* G4ExtDEDXTable::GetPhysicsVector(G4int atomicNumberIon, // Atomic number of ion
126 const G4String& matIdentifier // Name or chemical formula of material
127 )
128 {
129 G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
130
131 auto iter = dedxMapMaterials.find(key);
132
133 return (iter != dedxMapMaterials.end()) ? iter->second : nullptr;
134 }
135
136 // #########################################################################
137
138 G4double G4ExtDEDXTable::GetDEDX(G4double kinEnergyPerNucleon, // Kinetic energy per nucleon
139 G4int atomicNumberIon, // Atomic number of ion
140 G4int atomicNumberElem // Atomic number of elemental material
141 )
142 {
143 G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
144
145 auto iter = dedxMapElements.find(key);
146
147 return (iter != dedxMapElements.end()) ? (iter->second)->Value(kinEnergyPerNucleon) : 0.0;
148 }
149
150 // #########################################################################
151
152 G4double G4ExtDEDXTable::GetDEDX(G4double kinEnergyPerNucleon, // Kinetic energy per nucleon
153 G4int atomicNumberIon, // Atomic number of ion
154 const G4String& matIdentifier // Name or chemical formula of material
155 )
156 {
157 G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
158
159 auto iter = dedxMapMaterials.find(key);
160
161 return (iter != dedxMapMaterials.end()) ? (iter->second)->Value(kinEnergyPerNucleon) : 0.0;
162 }
163
164 // #########################################################################
165
166 G4bool G4ExtDEDXTable::AddPhysicsVector(G4PhysicsVector* physicsVector, // Physics vector
167 G4int atomicNumberIon, // Atomic number of ion
168 const G4String& matIdentifier, // Name of elemental material
169 G4int atomicNumberElem // Atomic number of elemental material
170 )
171 {
172 if (physicsVector == nullptr) {
173 G4Exception("G4ExtDEDXTable::AddPhysicsVector() for material", "mat037", FatalException,
174 "Pointer to vector is null-pointer.");
175 return false;
176 }
177
178 if (matIdentifier.empty()) {
179 G4Exception("G4ExtDEDXTable::AddPhysicsVector() for material", "mat038", FatalException,
180 "Invalid name of the material.");
181 return false;
182 }
183
184 if (atomicNumberIon <= 2) {
185 G4Exception("G4ExtDEDXTable::AddPhysicsVector() for material", "mat039", FatalException,
186 "Illegal atomic number.");
187 return false;
188 }
189
190 if (atomicNumberElem > 0) {
191 G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
192
193 if (dedxMapElements.count(key) == 1) {
194 G4Exception("G4ExtDEDXTable::AddPhysicsVector() for material", "mat037", FatalException,
195 "Vector already exist, remove it before replacing.");
196 return false;
197 }
198
199 dedxMapElements[key] = physicsVector;
200 }
201
202 G4IonDEDXKeyMat mkey = std::make_pair(atomicNumberIon, matIdentifier);
203
204 if (dedxMapMaterials.count(mkey) == 1) {
205 G4Exception("G4ExtDEDXTable::AddPhysicsVector() for material", "mat037", FatalException,
206 "Vector already exist, remove it before replacing.");
207 return false;
208 }
209
210 dedxMapMaterials[mkey] = physicsVector;
211
212 return true;
213 }
214
215 // #########################################################################
216
217 G4bool G4ExtDEDXTable::RemovePhysicsVector(G4int atomicNumberIon, // Atomic number of ion
218 const G4String& matIdentifier // Name or chemical formula of material
219 )
220 {
221 G4PhysicsVector* physicsVector = nullptr;
222
223 // Deleting key of physics vector from material map
224 G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
225
226 auto iter = dedxMapMaterials.find(key);
227
228 if (iter == dedxMapMaterials.end()) {
229 G4Exception("G4ExtDEDXTable::RemovePhysicsVector() for material", "mat037", FatalException,
230 "Pointer to vector is null-pointer.");
231 return false;
232 }
233
234 physicsVector = (*iter).second;
235 dedxMapMaterials.erase(key);
236
237 // Deleting key of physics vector from elemental material map (if it exists)
238 G4IonDEDXMapElem::iterator it;
239
240 for (it = dedxMapElements.begin(); it != dedxMapElements.end(); ++it) {
241 if ((*it).second == physicsVector) {
242 dedxMapElements.erase(it);
243 break;
244 }
245 }
246
247 // Deleting physics vector
248 delete physicsVector;
249
250 return true;
251 }
252
253 // #########################################################################
254
255 G4bool G4ExtDEDXTable::StorePhysicsTable(const G4String& fileName // File name
256 )
257 {
258 G4bool success = true;
259
260 std::ofstream ofilestream;
261
262 ofilestream.open(fileName, std::ios::out);
263
264 if (! ofilestream) {
265 G4ExceptionDescription ed;
266 ed << "Cannot open file " << fileName;
267 G4Exception("G4IonStoppingData::StorePhysicsTable()", "mat030", FatalException, ed);
268 success = false;
269 }
270 else {
271 size_t nmbMatTables = dedxMapMaterials.size();
272
273 ofilestream << nmbMatTables << G4endl << G4endl;
274
275 auto iterMat = dedxMapMaterials.begin();
276 auto iterMat_end = dedxMapMaterials.end();
277
278 for (; iterMat != iterMat_end; iterMat++) {
279 G4IonDEDXKeyMat key = iterMat->first;
280 G4PhysicsVector* physicsVector = iterMat->second;
281
282 G4int atomicNumberIon = key.first;
283 G4String matIdentifier = key.second;
284
285 G4int atomicNumberElem = FindAtomicNumberElement(physicsVector);
286
287 if (physicsVector != nullptr) {
288 ofilestream << atomicNumberIon << " " << matIdentifier;
289
290 if (atomicNumberElem > 0) {
291 ofilestream << " " << atomicNumberElem;
292 }
293
294 ofilestream << " # <Atomic number ion> <Material name> ";
295
296 if (atomicNumberElem > 0) {
297 ofilestream << "<Atomic number element>";
298 }
299
300 ofilestream << G4endl << physicsVector->GetType() << G4endl;
301
302 physicsVector->Store(ofilestream, true);
303
304 ofilestream << G4endl;
305 }
306 else {
307 G4Exception("G4IonStoppingData::StorePhysicsTable()", "mat030", FatalException,
308 "Cannot store vector.");
309 }
310 }
311 }
312
313 ofilestream.close();
314
315 return success;
316 }
317
318 // #########################################################################
319
320 G4bool G4ExtDEDXTable::RetrievePhysicsTable(const G4String& fileName)
321 {
322 std::ifstream ifilestream;
323 ifilestream.open(fileName, std::ios::in | std::ios::binary);
324 if (! ifilestream) {
325 G4ExceptionDescription ed;
326 ed << "Cannot open file " << fileName;
327 G4Exception("G4IonStoppingData::RetrievePhysicsTable()", "mat030", FatalException, ed);
328 return false;
329 }
330
331 // std::string::size_type nmbVectors;
332 G4int nmbVectors = 0;
333 ifilestream >> nmbVectors;
334 if (ifilestream.fail() || nmbVectors <= 0) {
335 G4cout << "G4ExtDEDXTable::RetrievePhysicsTable() "
336 << " File content of " << fileName << " ill-formated."
337 << " Nvectors= " << nmbVectors << G4endl;
338 ifilestream.close();
339 return false;
340 }
341
342 for (G4int i = 0; i < nmbVectors; ++i) {
343 G4String line = "";
344 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
345 while (line.empty()) {
346 getline(ifilestream, line);
347 if (ifilestream.fail()) {
348 G4cout << "G4ExtDEDXTable::RetrievePhysicsTable() "
349 << " File content of " << fileName << " ill-formated." << G4endl;
350 ifilestream.close();
351 return false;
352 }
353
354 std::string::size_type pos = line.find_first_of('#');
355 if (pos != std::string::npos && pos > 0) {
356 line = line.substr(0, pos);
357 }
358 }
359
360 std::istringstream headerstream(line);
361
362 std::string::size_type atomicNumberIon;
363 headerstream >> atomicNumberIon;
364
365 G4String materialName;
366 headerstream >> materialName;
367
368 if (headerstream.fail() || std::string::npos == atomicNumberIon) {
369 G4cout << "G4ExtDEDXTable::RetrievePhysicsTable() "
370 << " File content of " << fileName << " ill-formated "
371 << " (vector header)." << G4endl;
372 ifilestream.close();
373 return false;
374 }
375
376 std::string::size_type atomicNumberMat;
377 headerstream >> atomicNumberMat;
378
379 if (headerstream.eof() || std::string::npos == atomicNumberMat) {
380 atomicNumberMat = 0;
381 }
382
383 G4int vectorType;
384 ifilestream >> vectorType;
385
386 G4PhysicsVector* physicsVector = CreatePhysicsVector(vectorType);
387
388 if (physicsVector == nullptr) {
389 G4cout << "G4ExtDEDXTable::RetrievePhysicsTable "
390 << " illegal physics Vector type " << vectorType << " in " << fileName << G4endl;
391 ifilestream.close();
392 return false;
393 }
394
395 if (! physicsVector->Retrieve(ifilestream, true)) {
396 G4cout << "G4ExtDEDXTable::RetrievePhysicsTable() "
397 << " File content of " << fileName << " ill-formated." << G4endl;
398 ifilestream.close();
399 return false;
400 }
401 physicsVector->FillSecondDerivatives();
402
403 // Retrieved vector is added to material store
404 if (! AddPhysicsVector(
405 physicsVector, (G4int)atomicNumberIon, materialName, (G4int)atomicNumberMat))
406 {
407 delete physicsVector;
408 ifilestream.close();
409 return false;
410 }
411 }
412
413 ifilestream.close();
414
415 return true;
416 }
417
418 // #########################################################################
419
420 G4PhysicsVector* G4ExtDEDXTable::CreatePhysicsVector(G4int vectorType)
421 {
422 G4PhysicsVector* physicsVector = nullptr;
423
424 switch (vectorType) {
425 case T_G4PhysicsLinearVector:
426 physicsVector = new G4PhysicsLinearVector(true);
427 break;
428
429 case T_G4PhysicsLogVector:
430 physicsVector = new G4PhysicsLogVector(true);
431 break;
432
433 case T_G4PhysicsFreeVector:
434 physicsVector = new G4PhysicsFreeVector(true);
435 break;
436
437 default:
438 break;
439 }
440 return physicsVector;
441 }
442
443 // #########################################################################
444
445 G4int G4ExtDEDXTable::FindAtomicNumberElement(G4PhysicsVector* physicsVector)
446 {
447 G4int atomicNumber = 0;
448
449 auto iter = dedxMapElements.begin();
450 auto iter_end = dedxMapElements.end();
451
452 for (; iter != iter_end; ++iter) {
453 if ((*iter).second == physicsVector) {
454 G4IonDEDXKeyElem key = (*iter).first;
455 atomicNumber = key.second;
456 }
457 }
458
459 return atomicNumber;
460 }
461
462 // #########################################################################
463
464 void G4ExtDEDXTable::ClearTable()
465 {
466 auto iterMat = dedxMapMaterials.begin();
467 auto iterMat_end = dedxMapMaterials.end();
468
469 for (; iterMat != iterMat_end; ++iterMat) {
470 G4PhysicsVector* vec = iterMat->second;
471
472 delete vec;
473 }
474
475 dedxMapElements.clear();
476 dedxMapMaterials.clear();
477 }
478
479 // #########################################################################
480
481 void G4ExtDEDXTable::DumpMap()
482 {
483 auto iterMat = dedxMapMaterials.begin();
484 auto iterMat_end = dedxMapMaterials.end();
485
486 G4cout << std::setw(15) << std::right << "Atomic nmb ion" << std::setw(25) << std::right
487 << "Material name" << std::setw(25) << std::right << "Atomic nmb material" << G4endl;
488
489 for (; iterMat != iterMat_end; ++iterMat) {
490 G4IonDEDXKeyMat key = iterMat->first;
491 G4PhysicsVector* physicsVector = iterMat->second;
492
493 G4int atomicNumberIon = key.first;
494 G4String matIdentifier = key.second;
495
496 G4int atomicNumberElem = FindAtomicNumberElement(physicsVector);
497
498 if (physicsVector != nullptr) {
499 G4cout << std::setw(15) << std::right << atomicNumberIon << std::setw(25) << std::right
500 << matIdentifier << std::setw(25) << std::right;
501
502 if (atomicNumberElem > 0) {
503 G4cout << atomicNumberElem;
504 }
505 else {
506 G4cout << "N/A";
507 }
508
509 G4cout << G4endl;
510 }
511 }
512 }
513
514 // #########################################################################
515