| 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: G4IonStoppingData
30 //
31 // Base class: G4VIonDEDXTable
32 //
33 // Author: Anton Lechner (Anton.Lechner@cern.ch)
34 //
35 // First implementation: 03. 11. 2009
36 //
37 // Modifications:
38 // 25.10.2010 V.Ivanchenko fixed warnings reported by the Coverity tool
39 // 25.10.2011: new scheme for G4Exception (mma)
40 //
41 //
42 // Class description: Class which can read ion stopping power data from
43 // $G4LEDATA/ion_stopping_data
44 //
45 // Comments:
46 //
47 // ===========================================================================
48
49 #include "G4IonStoppingData.hh"
50
51 #include "G4PhysicalConstants.hh"
52 #include "G4PhysicsFreeVector.hh"
53 #include "G4PhysicsVector.hh"
54 #include "G4SystemOfUnits.hh"
55
56 #include <fstream>
57 #include <iomanip>
58 #include <sstream>
59 #include <utility>
60
61 // #########################################################################
62
63 G4IonStoppingData::G4IonStoppingData(const G4String& dir, G4bool val) : subDir(dir), fICRU90(val) {}
64
65 // #########################################################################
66
67 G4IonStoppingData::~G4IonStoppingData() { ClearTable(); }
68
69 // #########################################################################
70
71 G4bool G4IonStoppingData::IsApplicable(G4int atomicNumberIon, // Atomic number of ion
72 G4int atomicNumberElem // Atomic number of elemental material
73 )
74 {
75 G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
76
77 auto iter = dedxMapElements.find(key);
78
79 return iter != dedxMapElements.end();
80 }
81
82 // #########################################################################
83
84 G4bool G4IonStoppingData::IsApplicable(G4int atomicNumberIon, // Atomic number of ion
85 const G4String& matIdentifier // Name or chemical formula of material
86 )
87 {
88 G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
89
90 auto iter = dedxMapMaterials.find(key);
91
92 return iter != dedxMapMaterials.end();
93 }
94
95 // #########################################################################
96
97 G4PhysicsVector* G4IonStoppingData::GetPhysicsVector(G4int atomicNumberIon, // Atomic number of ion
98 G4int atomicNumberElem // Atomic number of elemental material
99 )
100 {
101 G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
102
103 auto iter = dedxMapElements.find(key);
104
105 return (iter != dedxMapElements.end()) ? iter->second : nullptr;
106 }
107
108 // #########################################################################
109
110 G4PhysicsVector* G4IonStoppingData::GetPhysicsVector(G4int atomicNumberIon, // Atomic number of ion
111 const G4String& matIdentifier // Name or chemical formula of material
112 )
113 {
114 G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
115
116 auto iter = dedxMapMaterials.find(key);
117
118 return (iter != dedxMapMaterials.end()) ? iter->second : nullptr;
119 }
120
121 // #########################################################################
122
123 G4double G4IonStoppingData::GetDEDX(G4double kinEnergyPerNucleon, // Kinetic energy per nucleon
124 G4int atomicNumberIon, // Atomic number of ion
125 G4int atomicNumberElem // Atomic number of elemental material
126 )
127 {
128 G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
129
130 auto iter = dedxMapElements.find(key);
131
132 return (iter != dedxMapElements.end()) ? (iter->second)->Value(kinEnergyPerNucleon) : 0.0;
133 }
134
135 // #########################################################################
136
137 G4double G4IonStoppingData::GetDEDX(G4double kinEnergyPerNucleon, // Kinetic energy per nucleon
138 G4int atomicNumberIon, // Atomic number of ion
139 const G4String& matIdentifier // Name or chemical formula of material
140 )
141 {
142 G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
143
144 auto iter = dedxMapMaterials.find(key);
145
146 return (iter != dedxMapMaterials.end()) ? (iter->second)->Value(kinEnergyPerNucleon) : 0.0;
147 }
148
149 // #########################################################################
150
151 G4bool G4IonStoppingData::AddPhysicsVector(G4PhysicsVector* physicsVector, // Physics vector
152 G4int atomicNumberIon, // Atomic number of ion
153 const G4String& matIdentifier // Name of elemental material
154 )
155 {
156 if (physicsVector == nullptr) {
157 G4Exception("G4IonStoppingData::AddPhysicsVector() for material", "mat037", FatalException,
158 "Pointer to vector is null-pointer.");
159 return false;
160 }
161
162 if (matIdentifier.empty()) {
163 G4Exception("G4IonStoppingData::AddPhysicsVector() for material", "mat038", FatalException,
164 "Invalid name of the material.");
165 return false;
166 }
167
168 if (atomicNumberIon <= 0) {
169 G4Exception("G4IonStoppingData::AddPhysicsVector() for material", "mat039", FatalException,
170 "Illegal atomic number.");
171 return false;
172 }
173
174 G4IonDEDXKeyMat mkey = std::make_pair(atomicNumberIon, matIdentifier);
175
176 if (dedxMapMaterials.count(mkey) == 1) {
177 G4ExceptionDescription ed;
178 ed << "Vector with Z1 = " << atomicNumberIon << ", mat = " << matIdentifier
179 << "already exists. Remove first before replacing.";
180 G4Exception("G4IonStoppingData::AddPhysicsVector() for material", "mat040", FatalException, ed);
181 return false;
182 }
183
184 dedxMapMaterials[mkey] = physicsVector;
185
186 return true;
187 }
188
189 // #########################################################################
190
191 G4bool G4IonStoppingData::AddPhysicsVector(G4PhysicsVector* physicsVector, // Physics vector
192 G4int atomicNumberIon, // Atomic number of ion
193 G4int atomicNumberElem // Atomic number of elemental material
194 )
195 {
196 if (physicsVector == nullptr) {
197 G4Exception("G4IonStoppingData::AddPhysicsVector() for element", "mat037", FatalException,
198 "Pointer to vector is null-pointer.");
199 return false;
200 }
201
202 if (atomicNumberIon <= 0) {
203 G4Exception("G4IonStoppingData::AddPhysicsVector() for element", "mat038", FatalException,
204 "Invalid ion number.");
205 return false;
206 }
207
208 if (atomicNumberElem <= 0) {
209 G4Exception("G4IonStoppingData::AddPhysicsVector() for element", "mat039", FatalException,
210 "Illegal atomic number.");
211 return false;
212 }
213
214 G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
215
216 if (dedxMapElements.count(key) == 1) {
217 G4ExceptionDescription ed;
218 ed << "Vector with Z1 = " << atomicNumberIon << ", Z= " << atomicNumberElem
219 << "already exists. Remove first before replacing.";
220 G4Exception("G4IonStoppingData::AddPhysicsVector() for element", "mat040", FatalException, ed);
221 return false;
222 }
223
224 dedxMapElements[key] = physicsVector;
225
226 return true;
227 }
228
229 // #########################################################################
230
231 G4bool G4IonStoppingData::RemovePhysicsVector(G4int atomicNumberIon, // Atomic number of ion
232 const G4String& matIdentifier // Name or chemical formula of material
233 )
234 {
235 G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
236
237 auto iter = dedxMapMaterials.find(key);
238
239 if (iter == dedxMapMaterials.end()) {
240 G4Exception("G4IonStoppingData::RemovePhysicsVector() for material", "mat038", FatalException,
241 "Invalid name of the material.");
242 return false;
243 }
244
245 G4PhysicsVector* physicsVector = (*iter).second;
246
247 // Deleting key of physics vector from material map
248 dedxMapMaterials.erase(key);
249
250 // Deleting physics vector
251 delete physicsVector;
252
253 return true;
254 }
255
256 // #########################################################################
257
258 G4bool G4IonStoppingData::RemovePhysicsVector(G4int atomicNumberIon, // Atomic number of ion
259 G4int atomicNumberElem // Atomic number of elemental material
260 )
261 {
262 G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
263
264 auto iter = dedxMapElements.find(key);
265
266 if (iter == dedxMapElements.end()) {
267 G4Exception("G4IonStoppingData::RemovePhysicsVector() for element", "mat038", FatalException,
268 "Invalid element.");
269 return false;
270 }
271
272 G4PhysicsVector* physicsVector = (*iter).second;
273
274 // Deleting key of physics vector from material map
275 dedxMapElements.erase(key);
276
277 // Deleting physics vector
278 delete physicsVector;
279
280 return true;
281 }
282
283 // #########################################################################
284
285 G4bool G4IonStoppingData::BuildPhysicsVector(G4int atomicNumberIon, // Atomic number of ion
286 const G4String& matname // Name of material
287 )
288 {
289 if (IsApplicable(atomicNumberIon, matname)) {
290 return true;
291 }
292
293 const char* path = G4FindDataDir("G4LEDATA");
294 if (path == nullptr) {
295 G4Exception("G4IonStoppingData::BuildPhysicsVector()", "mat521", FatalException,
296 "G4LEDATA environment variable not set");
297 return false;
298 }
299
300 std::ostringstream file;
301 G4String ww =
302 (fICRU90 && (matname == "G4_WATER" || matname == "G4_AIR" || matname == "G4_GRAPHITE")) ? "90"
303 : "73";
304
305 file << path << "/" << subDir << ww << "/z" << atomicNumberIon << "_" << matname << ".dat";
306 G4String fileName = G4String(file.str().c_str());
307
308 std::ifstream ifilestream(fileName);
309
310 if (! ifilestream.is_open()) {
311 return false;
312 }
313
314 auto* physicsVector = new G4PhysicsFreeVector(true);
315
316 if (! physicsVector->Retrieve(ifilestream, true)) {
317 ifilestream.close();
318 return false;
319 }
320
321 physicsVector->ScaleVector(MeV, MeV * cm2 / (0.001 * g));
322 physicsVector->FillSecondDerivatives();
323
324 // Retrieved vector is added to material store
325 if (! AddPhysicsVector(physicsVector, atomicNumberIon, matname)) {
326 delete physicsVector;
327 ifilestream.close();
328 return false;
329 }
330
331 ifilestream.close();
332 return true;
333 }
334
335 // #########################################################################
336
337 G4bool G4IonStoppingData::BuildPhysicsVector(G4int ZIon, // Atomic number of ion
338 G4int ZElem // Atomic number of elemental material
339 )
340 {
341 if (IsApplicable(ZIon, ZElem)) {
342 return true;
343 }
344
345 const char* path = G4FindDataDir("G4LEDATA");
346 if (path == nullptr) {
347 G4Exception("G4IonStoppingData::BuildPhysicsVector()", "mat522", FatalException,
348 "G4LEDATA environment variable not set");
349 return false;
350 }
351 std::ostringstream file;
352 G4String ww =
353 (fICRU90 && ZIon <= 18 && (ZElem == 1 || ZElem == 6 || ZElem == 7 || ZElem == 8)) ? "90" : "73";
354
355 file << path << "/" << subDir << ww << "/z" << ZIon << "_" << ZElem << ".dat";
356
357 G4String fileName = G4String(file.str().c_str());
358 std::ifstream ifilestream(fileName);
359
360 if (! ifilestream.is_open()) {
361 return false;
362 }
363 auto* physicsVector = new G4PhysicsFreeVector(true);
364
365 if (! physicsVector->Retrieve(ifilestream, true)) {
366 ifilestream.close();
367 return false;
368 }
369
370 physicsVector->ScaleVector(MeV, MeV * cm2 / (0.001 * g));
371 physicsVector->FillSecondDerivatives();
372
373 // Retrieved vector is added to material store
374 if (! AddPhysicsVector(physicsVector, ZIon, ZElem)) {
375 delete physicsVector;
376 ifilestream.close();
377 return false;
378 }
379
380 ifilestream.close();
381 return true;
382 }
383
384 // #########################################################################
385
386 void G4IonStoppingData::ClearTable()
387 {
388 auto iterMat = dedxMapMaterials.begin();
389 auto iterMat_end = dedxMapMaterials.end();
390
391 for (; iterMat != iterMat_end; iterMat++) {
392 G4PhysicsVector* vec = iterMat->second;
393
394 delete vec;
395 }
396
397 dedxMapMaterials.clear();
398
399 auto iterElem = dedxMapElements.begin();
400 auto iterElem_end = dedxMapElements.end();
401
402 for (; iterElem != iterElem_end; iterElem++) {
403 G4PhysicsVector* vec = iterElem->second;
404
405 delete vec;
406 }
407
408 dedxMapElements.clear();
409 }
410
411 // #########################################################################
412
413 void G4IonStoppingData::DumpMap()
414 {
415 auto iterMat = dedxMapMaterials.begin();
416 auto iterMat_end = dedxMapMaterials.end();
417
418 G4cout << std::setw(15) << std::right << "Atomic nmb ion" << std::setw(25) << std::right
419 << "Material name" << G4endl;
420
421 for (; iterMat != iterMat_end; iterMat++) {
422 G4IonDEDXKeyMat key = iterMat->first;
423 G4PhysicsVector* physicsVector = iterMat->second;
424
425 G4int atomicNumberIon = key.first;
426 G4String matIdentifier = key.second;
427
428 if (physicsVector != nullptr) {
429 G4cout << std::setw(15) << std::right << atomicNumberIon << std::setw(25) << std::right
430 << matIdentifier << G4endl;
431 }
432 }
433
434 auto iterElem = dedxMapElements.begin();
435 auto iterElem_end = dedxMapElements.end();
436
437 G4cout << std::setw(15) << std::right << "Atomic nmb ion" << std::setw(25) << std::right
438 << "Atomic nmb material" << G4endl;
439
440 for (; iterElem != iterElem_end; iterElem++) {
441 G4IonDEDXKeyElem key = iterElem->first;
442 G4PhysicsVector* physicsVector = iterElem->second;
443
444 G4int atomicNumberIon = key.first;
445 G4int atomicNumberElem = key.second;
446
447 if (physicsVector != nullptr) {
448 G4cout << std::setw(15) << std::right << atomicNumberIon << std::setw(25) << std::right
449 << atomicNumberElem << G4endl;
450 }
451 }
452 }
453
454 // #########################################################################
455