ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4eIonisationParameters.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4eIonisationParameters.cc
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 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
29 //
30 // History:
31 // -----------
32 // 31 Jul 2001 MGP Created, with dummy implementation
33 // 12.09.01 V.Ivanchenko Add param and interpolation of parameters
34 // 04.10.01 V.Ivanchenko Add BindingEnergy method
35 // 25.10.01 MGP Many bug fixes, mostly related to the
36 // management of pointers
37 // 29.11.01 V.Ivanchenko New parametrisation + Excitation
38 // 30.05.02 V.Ivanchenko Format and names of the data files were
39 // chenged to "ion-..."
40 // 17.02.04 V.Ivanchenko Increase buffer size
41 // 23.03.09 L.Pandola Updated warning message
42 // 03.12.10 V.Ivanchenko Fixed memory leak in LoadData
43 //
44 // -------------------------------------------------------------------
45 
47 #include "G4SystemOfUnits.hh"
48 #include "G4VEMDataSet.hh"
49 #include "G4ShellEMDataSet.hh"
50 #include "G4EMDataSet.hh"
51 #include "G4CompositeEMDataSet.hh"
52 #include "G4LinInterpolation.hh"
53 #include "G4LogLogInterpolation.hh"
56 #include "G4Material.hh"
57 #include "G4DataVector.hh"
58 #include <fstream>
59 #include <sstream>
60 
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,G4VEMDataSet*,std::less<G4int> >::iterator pos;
74 
75  for (pos = param.begin(); pos != param.end(); ++pos)
76  {
77  G4VEMDataSet* dataSet = (*pos).second;
78  delete dataSet;
79  }
80 
81  for (pos = excit.begin(); pos != excit.end(); ++pos)
82  {
83  G4VEMDataSet* 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,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
98 
99  pos = param.find(id);
100  if (pos!= param.end()) {
101  G4VEMDataSet* dataSet = (*pos).second;
102  G4int nShells = dataSet->NumberOfComponents();
103 
104  if(shellIndex < nShells) {
105  const G4VEMDataSet* 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,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
128 
129  pos = excit.find(Z);
130  if (pos!= excit.end()) {
131  G4VEMDataSet* 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("G4eIonisationParameters::LoadData",
157  "em1001",FatalException,"Unable to find MaterialTable");
158 
160 
161  for (G4int mLocal=0; mLocal<nMaterials; mLocal++) {
162 
163  const G4Material* material= (*materialTable)[mLocal];
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 = std::getenv("G4LEDATA");
177  if (!path)
178  {
179  G4Exception("G4BremsstrahlungParameters::LoadData",
180  "em0006",FatalException,"G4LEDATA environment variable not set");
181  return;
182  }
183 
184  G4String pathString(path);
185  G4String path2("/ioni/ion-sp-");
186  pathString += path2;
187 
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");
205  G4Exception("G4eIonisationParameters::LoadData",
206  "em0003",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<G4VEMDataSet*> p;
221  for (size_t k=0; k<length; k++)
222  {
224  G4VEMDataSet* composite = new G4CompositeEMDataSet(inter,1.,1.);
225  p.push_back(composite);
226  }
227 
228  G4int shell = 0;
229  std::vector<G4DataVector*> a;
230  a.resize(length);
231  G4DataVector e;
232  G4bool isReady = false;
233  do {
234  file >> energy >> sum;
235  //Check if energy is valid
236  if (energy < -2)
237  {
238  G4String excep("invalid data file");
239  G4Exception("G4eIonisationParameters::LoadData",
240  "em0005",FatalException,excep);
241  return;
242  }
243 
244  if (energy == -2) break;
245 
246  if (energy > -1) {
247  // new energy
248  if(!isReady) {
249  isReady = true;
250  e.clear();
251  for (size_t j=0; j<length; ++j) {
252  a[j] = new G4DataVector();
253  }
254  }
255  e.push_back(energy);
256  a[0]->push_back(sum);
257  for (size_t j=1; j<length; ++j) {
258  G4double qRead;
259  file >> qRead;
260  a[j]->push_back(qRead);
261  }
262 
263  } else {
264 
265  // End of set for a shell, fill the map
266  for (size_t k=0; k<length; k++) {
267 
268  G4VDataSetAlgorithm* interp;
269  if(0 == k) { interp = new G4LinLogLogInterpolation(); }
270  else { interp = new G4LogLogInterpolation(); }
271 
272  G4DataVector* eVector = new G4DataVector;
273  size_t eSize = e.size();
274  for (size_t sLocal=0; sLocal<eSize; sLocal++) {
275  eVector->push_back(e[sLocal]);
276  }
277  G4VEMDataSet* set = new G4EMDataSet(shell,eVector,a[k],interp,1.,1.);
278 
279  p[k]->AddComponent(set);
280  }
281 
282  // next shell
283  ++shell;
284  isReady = false;
285  }
286  } while (energy > -2);
287 
288  file.close();
289 
290  for (size_t kk=0; kk<length; kk++)
291  {
292  G4int id = Z*100 + kk;
293  param[id] = p[kk];
294  }
295  }
296 
297  G4String pathString_a(path);
298  G4String name_a = pathString_a + G4String("/ioni/ion-ex-av.dat");
299  std::ifstream file_a(name_a);
300  std::filebuf* lsdp_a = file_a.rdbuf();
301  G4String pathString_b(path);
302  G4String name_b = pathString_b + G4String("/ioni/ion-ex-sig.dat");
303  std::ifstream file_b(name_b);
304  std::filebuf* lsdp_b = file_b.rdbuf();
305 
306  if (! (lsdp_a->is_open()) ) {
307  G4String excep = G4String("cannot open file ")
308  + name_a;
309  G4Exception("G4eIonisationParameters::LoadData",
310  "em0003",FatalException,excep);
311  }
312  if (! (lsdp_b->is_open()) ) {
313  G4String excep = G4String("cannot open file ")
314  + name_b;
315  G4Exception("G4eIonisationParameters::LoadData",
316  "em0003",FatalException,excep);
317  }
318 
319  // The file is organized into two columns:
320  // 1st column is the energy
321  // 2nd column is the corresponding value
322  // The file terminates with the pattern: -1 -1
323  // -2 -2
324 
325  G4double ener, ener1, sig, sig1;
326  G4int z = 0;
327 
328  G4DataVector e;
329  e.clear();
330  G4DataVector d;
331  d.clear();
332 
333  do {
334  file_a >> ener >> sig;
335  file_b >> ener1 >> sig1;
336 
337  if(ener != ener1) {
338  G4cout << "G4eIonisationParameters: problem in excitation data "
339  << "ener= " << ener
340  << " ener1= " << ener1
341  << G4endl;
342  }
343 
344  // End of file
345  if (ener == -2) {
346  break;
347 
348  // End of next element
349  } else if (ener == -1) {
350 
351  z++;
352  G4double Z = (G4double)z;
353 
354  // fill map if Z is used
355  if (activeZ.contains(Z)) {
356 
358  G4DataVector* eVector = new G4DataVector;
359  G4DataVector* dVector = new G4DataVector;
360  size_t eSize = e.size();
361  for (size_t sLocal2=0; sLocal2<eSize; sLocal2++) {
362  eVector->push_back(e[sLocal2]);
363  dVector->push_back(d[sLocal2]);
364  }
365  G4VEMDataSet* set = new G4EMDataSet(z,eVector,dVector,inter,1.,1.);
366  excit[z] = set;
367  }
368  e.clear();
369  d.clear();
370 
371  } else {
372 
373  e.push_back(ener*MeV);
374  d.push_back(sig1*sig*barn*MeV);
375  }
376  } while (ener != -2);
377 
378  file_a.close();
379 
380 }
381 
382 
384 {
385  G4cout << G4endl;
386  G4cout << "===== G4eIonisationParameters =====" << G4endl;
387  G4cout << G4endl;
388 
389  size_t nZ = activeZ.size();
390  std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
391 
392  for (size_t i=0; i<nZ; i++) {
393  G4int Z = (G4int)activeZ[i];
394 
395  for (size_t j=0; j<length; j++) {
396 
397  G4int index = Z*100 + j;
398 
399  pos = param.find(index);
400  if (pos!= param.end()) {
401  G4VEMDataSet* dataSet = (*pos).second;
402  size_t nShells = dataSet->NumberOfComponents();
403 
404  for (size_t k=0; k<nShells; k++) {
405 
406  G4cout << "===== Z= " << Z << " shell= " << k
407  << " parameter[" << j << "] ====="
408  << G4endl;
409  const G4VEMDataSet* comp = dataSet->GetComponent(k);
410  comp->PrintData();
411  }
412  }
413  }
414  }
415  G4cout << "====================================" << G4endl;
416 }
417 
418