Geant4  10.01.p01
G4RDeIonisationParameters.cc
Go to the documentation of this file.
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 // $Id$
28 // GEANT4 tag $Name: geant4-09-01-ref-00 $
29 //
30 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
31 //
32 // History:
33 // -----------
34 // 31 Jul 2001 MGP Created, with dummy implementation
35 // 12.09.01 V.Ivanchenko Add param and interpolation of parameters
36 // 04.10.01 V.Ivanchenko Add BindingEnergy method
37 // 25.10.01 MGP Many bug fixes, mostly related to the
38 // management of pointers
39 // 29.11.01 V.Ivanchenko New parametrisation + Excitation
40 // 30.05.02 V.Ivanchenko Format and names of the data files were
41 // chenged to "ion-..."
42 // 17.02.04 V.Ivanchenko Increase buffer size
43 //
44 // -------------------------------------------------------------------
45 
46 #include <fstream>
47 #include <sstream>
48 
50 #include "G4RDVEMDataSet.hh"
51 #include "G4RDShellEMDataSet.hh"
52 #include "G4RDEMDataSet.hh"
54 #include "G4RDLinInterpolation.hh"
58 #include "G4SystemOfUnits.hh"
59 #include "G4Material.hh"
60 #include "G4DataVector.hh"
61 
63  : zMin(minZ), zMax(maxZ),
64  length(24)
65 {
66  LoadData();
67 }
68 
69 
71 {
72  // Reset the map of data sets: remove the data sets from the map
73  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::iterator pos;
74 
75  for (pos = param.begin(); pos != param.end(); ++pos)
76  {
77  G4RDVEMDataSet* dataSet = (*pos).second;
78  delete dataSet;
79  }
80 
81  for (pos = excit.begin(); pos != excit.end(); ++pos)
82  {
83  G4RDVEMDataSet* dataSet = (*pos).second;
84  delete dataSet;
85  }
86 
87  activeZ.clear();
88 }
89 
90 
92  G4int parameterIndex,
93  G4double e) const
94 {
95  G4double value = 0.;
96  G4int id = Z*100 + parameterIndex;
97  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
98 
99  pos = param.find(id);
100  if (pos!= param.end()) {
101  G4RDVEMDataSet* dataSet = (*pos).second;
102  G4int nShells = dataSet->NumberOfComponents();
103 
104  if(shellIndex < nShells) {
105  const G4RDVEMDataSet* component = dataSet->GetComponent(shellIndex);
106  const G4DataVector ener = component->GetEnergies(0);
107  G4double ee = std::max(ener.front(),std::min(ener.back(),e));
108  value = component->FindValue(ee);
109  } else {
110  G4cout << "WARNING: G4IonisationParameters::FindParameter "
111  << "has no parameters for shell= " << shellIndex
112  << "; Z= " << Z
113  << G4endl;
114  }
115  } else {
116  G4cout << "WARNING: G4IonisationParameters::Parameter "
117  << "did not find ID = "
118  << shellIndex << G4endl;
119  }
120 
121  return value;
122 }
123 
125 {
126  G4double value = 0.;
127  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
128 
129  pos = excit.find(Z);
130  if (pos!= excit.end()) {
131  G4RDVEMDataSet* dataSet = (*pos).second;
132 
133  const G4DataVector ener = dataSet->GetEnergies(0);
134  G4double ee = std::max(ener.front(),std::min(ener.back(),e));
135  value = dataSet->FindValue(ee);
136  } else {
137  G4cout << "WARNING: G4IonisationParameters::Excitation "
138  << "did not find ID = "
139  << Z << G4endl;
140  }
141 
142  return value;
143 }
144 
145 
147 {
148  // ---------------------------------------
149  // Please document what are the parameters
150  // ---------------------------------------
151 
152  // define active elements
153 
154  const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
155  if (materialTable == 0)
156  G4Exception("G4RDeIonisationParameters::LoadData()", "InvalidSetup",
157  FatalException, "No MaterialTable found!");
158 
160 
161  for (G4int m=0; m<nMaterials; m++) {
162 
163  const G4Material* material= (*materialTable)[m];
164  const G4ElementVector* elementVector = material->GetElementVector();
165  const size_t nElements = material->GetNumberOfElements();
166 
167  for (size_t iEl=0; iEl<nElements; iEl++) {
168  G4Element* element = (*elementVector)[iEl];
169  G4double Z = element->GetZ();
170  if (!(activeZ.contains(Z))) {
171  activeZ.push_back(Z);
172  }
173  }
174  }
175 
176  char* path = getenv("G4LEDATA");
177  if (!path)
178  {
179  G4String excep("G4LEDATA environment variable not set!");
180  G4Exception("G4RDeIonisationParameters::LoadData()", "InvalidSetup",
181  FatalException, excep);
182  }
183 
184  G4String pathString(path);
185  G4String path2("/ioni/ion-sp-");
186  pathString += path2;
187 
188  G4double energy, sum;
189 
190  size_t nZ = activeZ.size();
191 
192  for (size_t i=0; i<nZ; i++) {
193 
194  G4int Z = (G4int)activeZ[i];
195  std::ostringstream ost;
196  ost << pathString << Z << ".dat";
197  G4String name(ost.str());
198 
199  std::ifstream file(name);
200  std::filebuf* lsdp = file.rdbuf();
201 
202  if (! (lsdp->is_open()) ) {
203  G4String excep = G4String("Data file: ")
204  + name + G4String(" not found. The version 1.# of G4LEDATA should be used");
205  G4Exception("G4RDeIonisationParameters::LoadData()", "DataNotFound",
206  FatalException, excep);
207  }
208 
209  // The file is organized into:
210  // For each shell there are two lines:
211  // 1st line:
212  // 1st column is the energy of incident e-,
213  // 2d column is the parameter of screan term;
214  // 2d line:
215  // 3 energy (MeV) subdividing different approximation area of the spectrum
216  // 20 point on the spectrum
217  // The shell terminates with the pattern: -1 -1
218  // The file terminates with the pattern: -2 -2
219 
220  std::vector<G4RDVEMDataSet*> p;
221  for (size_t k=0; k<length; k++)
222  {
224  G4RDVEMDataSet* composite = new G4RDCompositeEMDataSet(inter,1.,1.);
225  p.push_back(composite);
226  }
227 
228  G4int shell = 0;
229  std::vector<G4DataVector*> a;
230  for (size_t j=0; j<length; j++)
231  {
232  G4DataVector* aa = new G4DataVector();
233  a.push_back(aa);
234  }
235  G4DataVector e;
236  e.clear();
237  do {
238  file >> energy >> sum;
239  if (energy == -2) break;
240 
241  if (energy > -1) {
242  e.push_back(energy);
243  a[0]->push_back(sum);
244  for (size_t j=0; j<length-1; j++) {
245  G4double qRead;
246  file >> qRead;
247  a[j + 1]->push_back(qRead);
248  }
249 
250  } else {
251 
252  // End of set for a shell, fill the map
253  for (size_t k=0; k<length; k++) {
254 
255  G4RDVDataSetAlgorithm* interp;
256  if(0 == k) interp = new G4RDLinLogLogInterpolation();
257  else interp = new G4RDLogLogInterpolation();
258 
259  G4DataVector* eVector = new G4DataVector;
260  size_t eSize = e.size();
261  for (size_t s=0; s<eSize; s++) {
262  eVector->push_back(e[s]);
263  }
264  G4RDVEMDataSet* set = new G4RDEMDataSet(shell,eVector,a[k],interp,1.,1.);
265 
266  p[k]->AddComponent(set);
267  }
268 
269  // clear vectors
270  for (size_t j2=0; j2<length; j2++) {
271  a[j2] = new G4DataVector();
272  }
273  shell++;
274  e.clear();
275  }
276  } while (energy > -2);
277 
278  file.close();
279 
280  for (size_t kk=0; kk<length; kk++)
281  {
282  G4int id = Z*100 + kk;
283  param[id] = p[kk];
284  }
285  }
286 
287  G4String pathString_a(path);
288  G4String name_a = pathString_a + G4String("/ioni/ion-ex-av.dat");
289  std::ifstream file_a(name_a);
290  std::filebuf* lsdp_a = file_a.rdbuf();
291  G4String pathString_b(path);
292  G4String name_b = pathString_b + G4String("/ioni/ion-ex-sig.dat");
293  std::ifstream file_b(name_b);
294  std::filebuf* lsdp_b = file_b.rdbuf();
295 
296  if (! (lsdp_a->is_open()) ) {
297  G4String excep = G4String("Cannot open file ")
298  + name_a;
299  G4Exception("G4RDeIonisationParameters::LoadData()", "CannotOpenFile",
300  FatalException, excep);
301  }
302  if (! (lsdp_b->is_open()) ) {
303  G4String excep = G4String("Cannot open file ")
304  + name_b;
305  G4Exception("G4RDeIonisationParameters::LoadData()", "CannotOpenFile",
306  FatalException, excep);
307  }
308 
309  // The file is organized into two columns:
310  // 1st column is the energy
311  // 2nd column is the corresponding value
312  // The file terminates with the pattern: -1 -1
313  // -2 -2
314 
315  G4double ener, ener1, sig, sig1;
316  G4int z = 0;
317 
318  G4DataVector e;
319  e.clear();
320  G4DataVector d;
321  d.clear();
322 
323  do {
324  file_a >> ener >> sig;
325  file_b >> ener1 >> sig1;
326 
327  if(ener != ener1) {
328  G4cout << "G4RDeIonisationParameters: problem in excitation data "
329  << "ener= " << ener
330  << " ener1= " << ener1
331  << G4endl;
332  }
333 
334  // End of file
335  if (ener == -2) {
336  break;
337 
338  // End of next element
339  } else if (ener == -1) {
340 
341  z++;
342  G4double Z = (G4double)z;
343 
344  // fill map if Z is used
345  if (activeZ.contains(Z)) {
346 
348  G4DataVector* eVector = new G4DataVector;
349  G4DataVector* dVector = new G4DataVector;
350  size_t eSize = e.size();
351  for (size_t s=0; s<eSize; s++) {
352  eVector->push_back(e[s]);
353  dVector->push_back(d[s]);
354  }
355  G4RDVEMDataSet* set = new G4RDEMDataSet(z,eVector,dVector,inter,1.,1.);
356  excit[z] = set;
357  }
358  e.clear();
359  d.clear();
360 
361  } else {
362 
363  e.push_back(ener*MeV);
364  d.push_back(sig1*sig*barn*MeV);
365  }
366  } while (ener != -2);
367 
368  file_a.close();
369 
370 }
371 
372 
374 {
375  G4cout << G4endl;
376  G4cout << "===== G4RDeIonisationParameters =====" << G4endl;
377  G4cout << G4endl;
378 
379  size_t nZ = activeZ.size();
380  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
381 
382  for (size_t i=0; i<nZ; i++) {
383  G4int Z = (G4int)activeZ[i];
384 
385  for (size_t j=0; j<length; j++) {
386 
387  G4int index = Z*100 + j;
388 
389  pos = param.find(index);
390  if (pos!= param.end()) {
391  G4RDVEMDataSet* dataSet = (*pos).second;
392  size_t nShells = dataSet->NumberOfComponents();
393 
394  for (size_t k=0; k<nShells; k++) {
395 
396  G4cout << "===== Z= " << Z << " shell= " << k
397  << " parameter[" << j << "] ====="
398  << G4endl;
399  const G4RDVEMDataSet* comp = dataSet->GetComponent(k);
400  comp->PrintData();
401  }
402  }
403  }
404  }
405  G4cout << "====================================" << G4endl;
406 }
407 
408 
static const double MeV
Definition: G4SIunits.hh:193
std::vector< G4Element * > G4ElementVector
G4double z
Definition: TRTMaterials.hh:39
G4String name
Definition: TRTMaterials.hh:40
G4double GetZ() const
Definition: G4Element.hh:131
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:588
std::vector< G4Material * > G4MaterialTable
virtual const G4RDVEMDataSet * GetComponent(G4int componentId) const =0
G4double Excitation(G4int Z, G4double e) const
G4double a
Definition: TRTMaterials.hh:39
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
virtual size_t NumberOfComponents(void) const =0
static const double s
Definition: G4SIunits.hh:150
G4GLOB_DLL std::ostream G4cout
std::map< G4int, G4RDVEMDataSet *, std::less< G4int > > excit
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
virtual void PrintData(void) const =0
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:595
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual const G4DataVector & GetEnergies(G4int componentId) const =0
std::map< G4int, G4RDVEMDataSet *, std::less< G4int > > param
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
#define G4endl
Definition: G4ios.hh:61
static const double m
Definition: G4SIunits.hh:110
G4double Parameter(G4int Z, G4int shellIndex, G4int parameterIndex, G4double e) const
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
static const double barn
Definition: G4SIunits.hh:95
double G4double
Definition: G4Types.hh:76
G4bool contains(const G4double &) const
static const G4double pos
G4RDeIonisationParameters(G4int minZ=1, G4int maxZ=99)