source: trunk/source/materials/src/G4IonStoppingData.cc @ 1197

Last change on this file since 1197 was 1197, checked in by garnier, 15 years ago

nvx fichiers dans CVS
File size: 15.6 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:                G4IonStoppingData
31//
32// Base class:           G4VIonDEDXTable
33//
34// Author:               Anton Lechner (Anton.Lechner@cern.ch)
35//
36// First implementation: 03. 11. 2009
37//
38// Modifications:
39//
40//
41// Class description: Class which can read ion stopping power data from
42//                    $G4LEDATA/ion_stopping_data
43//
44// Comments:
45//
46// ===========================================================================
47//
48
49#include "G4IonStoppingData.hh"
50#include "G4PhysicsVector.hh"
51#include "G4LPhysicsFreeVector.hh"
52#include <fstream>
53#include <sstream>
54#include <iomanip>
55
56
57// #########################################################################
58
59G4IonStoppingData::G4IonStoppingData(const G4String& leDirectory) :
60  subDir( leDirectory ) {
61
62}
63
64// #########################################################################
65
66G4IonStoppingData::~G4IonStoppingData() {
67
68  ClearTable();
69}
70
71// #########################################################################
72
73G4bool G4IonStoppingData::IsApplicable(
74         G4int atomicNumberIon,  // Atomic number of ion
75         G4int atomicNumberElem  // Atomic number of elemental material
76                                    ) {
77  G4bool isApplicable = true; 
78  G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
79
80  G4IonDEDXMapElem::iterator iter = dedxMapElements.find(key);
81
82  if(iter == dedxMapElements.end()) isApplicable = false; 
83
84  return isApplicable; 
85}
86
87// #########################################################################
88
89G4bool G4IonStoppingData::IsApplicable(
90         G4int atomicNumberIon,         // Atomic number of ion
91         const G4String& matIdentifier  // Name or chemical formula of material
92                                    ) {
93  G4bool isApplicable = true; 
94  G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
95
96  G4IonDEDXMapMat::iterator iter = dedxMapMaterials.find(key);
97
98  if(iter == dedxMapMaterials.end()) isApplicable = false; 
99
100  return isApplicable; 
101}
102
103// #########################################################################
104
105G4PhysicsVector* G4IonStoppingData::GetPhysicsVector(
106         G4int atomicNumberIon,        // Atomic number of ion
107         G4int atomicNumberElem        // Atomic number of elemental material
108                                    ) {
109
110  G4PhysicsVector* physVector = 0;
111
112  G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
113
114  G4IonDEDXMapElem::iterator iter = dedxMapElements.find(key);
115
116  if(iter != dedxMapElements.end()) physVector = iter -> second; 
117
118  return physVector; 
119}
120
121// #########################################################################
122
123G4PhysicsVector*  G4IonStoppingData::GetPhysicsVector(
124         G4int atomicNumberIon,        // Atomic number of ion
125         const G4String& matIdentifier // Name or chemical formula of material
126                                    ) {
127
128  G4PhysicsVector* physVector = 0;
129
130  G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
131
132  G4IonDEDXMapMat::iterator iter = dedxMapMaterials.find(key);
133
134  if(iter != dedxMapMaterials.end()) physVector = iter -> second; 
135
136  return physVector; 
137}
138
139// #########################################################################
140
141G4double G4IonStoppingData::GetDEDX(
142         G4double kinEnergyPerNucleon, // Kinetic energy per nucleon
143         G4int atomicNumberIon,        // Atomic number of ion
144         G4int atomicNumberElem        // Atomic number of elemental material
145                                  ) {
146  G4double dedx = 0;
147
148  G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
149
150  G4IonDEDXMapElem::iterator iter = dedxMapElements.find(key);
151
152  if( iter != dedxMapElements.end() ) {
153     G4PhysicsVector* physicsVector = iter -> second; 
154
155     G4bool b;
156     dedx = physicsVector -> GetValue( kinEnergyPerNucleon, b );   
157  }
158
159  return dedx; 
160}
161
162// #########################################################################
163
164G4double G4IonStoppingData::GetDEDX(
165         G4double kinEnergyPerNucleon, // Kinetic energy per nucleon
166         G4int atomicNumberIon,        // Atomic number of ion
167         const G4String& matIdentifier // Name or chemical formula of material
168                                  ) {
169  G4double dedx = 0;
170
171  G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
172
173  G4IonDEDXMapMat::iterator iter = dedxMapMaterials.find(key);
174
175  if(iter != dedxMapMaterials.end()) {
176     G4PhysicsVector* physicsVector = iter -> second; 
177
178     G4bool b;
179     dedx = physicsVector -> GetValue( kinEnergyPerNucleon, b );   
180  }
181
182  return dedx; 
183}
184
185// #########################################################################
186
187G4bool G4IonStoppingData::AddPhysicsVector(
188        G4PhysicsVector* physicsVector, // Physics vector
189        G4int atomicNumberIon,          // Atomic number of ion
190        const G4String& matIdentifier   // Name of elemental material
191                                      ) {
192
193  if(physicsVector == 0) {
194
195#ifdef G4VERBOSE
196     G4cerr << "G4IonStoppingData::AddPhysicsVector() Error: Pointer to vector"
197            << " is null-pointer."
198            << G4endl;
199#endif
200
201     return false;
202  }
203
204  if(matIdentifier.empty()) {
205
206#ifdef G4VERBOSE
207     G4cerr << "G4IonStoppingData::AddPhysicsVector() Error: "
208            << "Cannot add physics vector. Invalid name."
209            << G4endl;
210#endif
211
212     return false;
213  }
214
215  if(atomicNumberIon <= 0) {
216
217#ifdef G4VERBOSE
218     G4cerr << "G4IonStoppingData::AddPhysicsVector() Error: "
219            << "Cannot add physics vector. Illegal atomic number."
220            << G4endl;
221#endif
222
223     return false;
224  }
225
226  G4IonDEDXKeyMat mkey = std::make_pair(atomicNumberIon, matIdentifier);
227
228  if(dedxMapMaterials.count(mkey) == 1) {
229
230#ifdef G4VERBOSE
231     G4cerr << "G4IonStoppingData::AddPhysicsVector() Error: "
232            << "Vector with Z1 = " << atomicNumberIon << ", mat = " 
233            << matIdentifier
234            << "already exists. Remove first before replacing."
235            << G4endl;
236#endif
237
238     return false;
239  }
240
241  dedxMapMaterials[mkey] = physicsVector;
242
243  return true;
244}
245
246// #########################################################################
247
248G4bool G4IonStoppingData::AddPhysicsVector(
249        G4PhysicsVector* physicsVector, // Physics vector
250        G4int atomicNumberIon,          // Atomic number of ion
251        G4int atomicNumberElem          // Atomic number of elemental material
252                                      ) {
253
254  if(physicsVector == 0) {
255
256#ifdef G4VERBOSE
257     G4cerr << "G4IonStoppingData::AddPhysicsVector() Error: "
258            << "Pointer to vector is null-pointer."
259            << G4endl;
260#endif
261
262     return false;
263  }
264
265  if(atomicNumberIon <= 0) {
266
267#ifdef G4VERBOSE
268     G4cerr << "G4IonStoppingData::AddPhysicsVector() Error: "
269            << "Cannot add physics vector. Illegal atomic number."
270            << G4endl;
271#endif
272
273     return false;
274  }
275
276  if(atomicNumberElem <= 0) {
277
278#ifdef G4VERBOSE
279        G4cerr << "G4IonStoppingData::AddPhysicsVector() Error: "
280               << "Atomic number of element < 0."
281               << G4endl;
282#endif
283        return false;
284  }
285
286  G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
287
288  if(dedxMapElements.count(key) == 1) {
289
290#ifdef G4VERBOSE
291     G4cerr << "G4IonStoppingData::AddPhysicsVector() Error: "
292            << "Vector with Z1 = " << atomicNumberIon << ", Z2 = " 
293            << atomicNumberElem
294            << " already exists. Remove first before replacing."
295            << G4endl;
296#endif
297      return false;
298  }
299
300  dedxMapElements[key] = physicsVector;
301
302  return true;
303}
304
305// #########################################################################
306
307G4bool G4IonStoppingData::RemovePhysicsVector(
308        G4int atomicNumberIon,         // Atomic number of ion
309        const G4String& matIdentifier  // Name or chemical formula of material
310                                      ) {
311
312  G4IonDEDXKeyMat key = std::make_pair(atomicNumberIon, matIdentifier);
313
314  G4IonDEDXMapMat::iterator iter = dedxMapMaterials.find(key);
315
316  if(iter == dedxMapMaterials.end()) {
317
318#ifdef G4VERBOSE
319     G4cerr << "G4IonStoppingData::RemovePhysicsVector() Warning: "
320            << "Cannot remove physics vector. Vector not found."
321            << G4endl;
322#endif
323
324     return false;
325  }
326
327  G4PhysicsVector* physicsVector = (*iter).second;
328
329  // Deleting key of physics vector from material map
330  dedxMapMaterials.erase(key);
331
332  // Deleting physics vector
333  delete physicsVector;
334
335  return true;
336}
337
338// #########################################################################
339
340G4bool G4IonStoppingData::RemovePhysicsVector(
341        G4int atomicNumberIon,         // Atomic number of ion
342        G4int atomicNumberElem         // Atomic number of elemental material
343                                      ) {
344  G4IonDEDXKeyElem key = std::make_pair(atomicNumberIon, atomicNumberElem);
345
346  G4IonDEDXMapElem::iterator iter = dedxMapElements.find(key);
347
348  if(iter == dedxMapElements.end()) {
349
350#ifdef G4VERBOSE
351     G4cerr << "G4IonStoppingData::RemovePhysicsVector() Warning: "
352            << "Cannot remove physics vector. Vector not found."
353            << G4endl;
354#endif
355
356     return false;
357  }
358
359  G4PhysicsVector* physicsVector = (*iter).second;
360
361  // Deleting key of physics vector from material map
362  dedxMapElements.erase(key);
363
364  // Deleting physics vector
365  delete physicsVector;
366
367  return true;
368}
369
370// #########################################################################
371
372G4bool G4IonStoppingData::BuildPhysicsVector(
373        G4int atomicNumberIon,          // Atomic number of ion
374        const G4String& matIdentifier   // Name of material
375                                                     ) {
376
377  if( IsApplicable(atomicNumberIon, matIdentifier) ) return true;
378
379  char* path = getenv("G4LEDATA");
380  if ( !path ) G4Exception("G4LEDATA environment variable not set");
381 
382  std::ostringstream file;
383 
384  file << path << "/" << subDir << "/z" 
385       << atomicNumberIon << "_" << matIdentifier
386       << ".dat";
387                     
388  G4String fileName = G4String( file.str().c_str() );
389
390  std::ifstream ifilestream( fileName );
391
392  if ( !ifilestream.is_open() ) return false;
393 
394  G4LPhysicsFreeVector* physicsVector = new G4LPhysicsFreeVector(); 
395
396  if( !physicsVector -> Retrieve(ifilestream, true) ) {
397 
398     ifilestream.close();
399     return false;
400  }
401
402  physicsVector -> ScaleVector( MeV, MeV * cm2 /( 0.001 * g) ); 
403  physicsVector -> SetSpline( true );
404  physicsVector -> FillSecondDerivatives();
405
406  // Retrieved vector is added to material store
407  if( !AddPhysicsVector(physicsVector, atomicNumberIon, matIdentifier) ) {
408     delete physicsVector;
409     ifilestream.close();
410     return false;
411  }
412
413  ifilestream.close();
414  return true;
415}
416
417// #########################################################################
418
419G4bool G4IonStoppingData::BuildPhysicsVector(
420        G4int atomicNumberIon,          // Atomic number of ion
421        G4int atomicNumberElem          // Atomic number of elemental material
422                                                     ) {
423
424  if( IsApplicable(atomicNumberIon, atomicNumberElem) ) return true;
425
426  char* path = getenv("G4LEDATA");
427  if ( !path ) G4Exception("G4LEDATA environment variable not set");
428 
429  std::ostringstream file;
430 
431  file << path << "/" << subDir << "/z" 
432       << atomicNumberIon << "_" << atomicNumberElem
433       << ".dat";
434                     
435  G4String fileName = G4String( file.str().c_str() );
436
437  std::ifstream ifilestream( fileName );
438
439  if ( !ifilestream.is_open() ) return false;
440 
441  G4LPhysicsFreeVector* physicsVector = new G4LPhysicsFreeVector(); 
442
443  if( !physicsVector -> Retrieve(ifilestream, true) ) {
444 
445     ifilestream.close();
446     return false;
447  }
448
449  physicsVector -> ScaleVector( MeV, MeV * cm2 /( 0.001 * g) ); 
450  physicsVector -> SetSpline( true );
451  physicsVector -> FillSecondDerivatives();
452
453  // Retrieved vector is added to material store
454  if( !AddPhysicsVector(physicsVector, atomicNumberIon, atomicNumberElem) ) {
455     delete physicsVector;
456     ifilestream.close();
457     return false;
458  }
459
460  ifilestream.close();
461  return true;
462}
463
464// #########################################################################
465
466void G4IonStoppingData::ClearTable() {
467
468  G4IonDEDXMapMat::iterator iterMat = dedxMapMaterials.begin();
469  G4IonDEDXMapMat::iterator iterMat_end = dedxMapMaterials.end();
470
471  for(;iterMat != iterMat_end; iterMat++) { 
472
473    G4PhysicsVector* vec = iterMat -> second;
474
475    if(vec != 0) delete vec;
476  }
477
478  dedxMapMaterials.clear();
479
480  G4IonDEDXMapElem::iterator iterElem = dedxMapElements.begin();
481  G4IonDEDXMapElem::iterator iterElem_end = dedxMapElements.end();
482
483  for(;iterElem != iterElem_end; iterElem++) { 
484
485    G4PhysicsVector* vec = iterElem -> second;
486
487    if(vec != 0) delete vec;
488  }
489
490  dedxMapElements.clear();
491}
492
493// #########################################################################
494
495void G4IonStoppingData::DumpMap() {
496
497  G4IonDEDXMapMat::iterator iterMat = dedxMapMaterials.begin();
498  G4IonDEDXMapMat::iterator iterMat_end = dedxMapMaterials.end();
499
500  G4cout << std::setw(15) << std::right
501         << "Atomic nmb ion"
502         << std::setw(25) << std::right
503         << "Material name"
504         << G4endl;
505
506  for(;iterMat != iterMat_end; iterMat++) {
507      G4IonDEDXKeyMat key = iterMat -> first;
508      G4PhysicsVector* physicsVector = iterMat -> second; 
509
510      G4int atomicNumberIon = key.first;
511      G4String matIdentifier = key.second;
512
513      if(physicsVector != 0) {
514         G4cout << std::setw(15) << std::right
515                << atomicNumberIon
516                << std::setw(25) << std::right
517                << matIdentifier
518                << G4endl;
519      }
520  }
521
522  G4IonDEDXMapElem::iterator iterElem = dedxMapElements.begin();
523  G4IonDEDXMapElem::iterator iterElem_end = dedxMapElements.end();
524
525  G4cout << std::setw(15) << std::right
526         << "Atomic nmb ion"
527         << std::setw(25) << std::right
528         << "Atomic nmb material"
529         << G4endl;
530
531  for(;iterElem != iterElem_end; iterElem++) { 
532      G4IonDEDXKeyElem key = iterElem -> first;
533      G4PhysicsVector* physicsVector = iterElem -> second; 
534
535      G4int atomicNumberIon = key.first;
536      G4int atomicNumberElem = key.second;
537
538      if(physicsVector != 0) {
539         G4cout << std::setw(15) << std::right
540                << atomicNumberIon
541                << std::setw(25) << std::right
542                << atomicNumberElem
543                << G4endl;
544      }
545  }
546
547}
548
549// #########################################################################
550
Note: See TracBrowser for help on using the repository browser.