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