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