| Geant4 Cross Reference |
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
46 #include "G4eIonisationParameters.hh"
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"
54 #include "G4LinLogLogInterpolation.hh"
55 #include "G4SemiLogInterpolation.hh"
56 #include "G4Material.hh"
57 #include "G4DataVector.hh"
58 #include <fstream>
59 #include <sstream>
60
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
62
63 G4eIonisationParameters:: G4eIonisationParameters()
64 {
65 LoadData();
66 }
67
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69
70 G4eIonisationParameters::~G4eIonisationParameters()
71 {
72 // Reset the map of data sets: remove the data sets from the map
73 for (auto& pos : param)
74 {
75 G4VEMDataSet* dataSet = pos.second;
76 delete dataSet;
77 dataSet = nullptr;
78 }
79 for (auto& pos : excit)
80 {
81 G4VEMDataSet* dataSet = pos.second;
82 delete dataSet;
83 dataSet = nullptr;
84 }
85 activeZ.clear();
86 }
87
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89
90 G4double G4eIonisationParameters::Parameter(G4int Z, G4int shellIndex,
91 G4int parameterIndex,
92 G4double e) const
93 {
94 G4double value = 0.;
95 G4int id = Z*100 + parameterIndex;
96
97 auto pos = param.find(id);
98 if (pos!= param.end()) {
99 G4VEMDataSet* dataSet = (*pos).second;
100 G4int nShells = (G4int)dataSet->NumberOfComponents();
101
102 if(shellIndex < nShells) {
103 const G4VEMDataSet* 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
122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123
124 G4double G4eIonisationParameters::Excitation(G4int Z, G4double e) const
125 {
126 G4double value = 0.;
127 auto pos = excit.find(Z);
128 if (pos!= excit.end()) {
129 G4VEMDataSet* 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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
144
145 void G4eIonisationParameters::LoadData()
146 {
147 // ---------------------------------------
148 // Please document what are the parameters
149 // ---------------------------------------
150
151 // define active elements
152
153 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
154 if (materialTable == nullptr)
155 G4Exception("G4eIonisationParameters::LoadData",
156 "em1001",FatalException,"Unable to find MaterialTable");
157
158 std::size_t nMaterials = G4Material::GetNumberOfMaterials();
159
160 for (std::size_t mLocal=0; mLocal<nMaterials; ++mLocal) {
161
162 const G4Material* material= (*materialTable)[mLocal];
163 const G4ElementVector* elementVector = material->GetElementVector();
164 const std::size_t nElements = material->GetNumberOfElements();
165
166 for (std::size_t iEl=0; iEl<nElements; ++iEl) {
167 G4double Z = (*elementVector)[iEl]->GetZ();
168 if (!(activeZ.contains(Z))) {
169 activeZ.push_back(Z);
170 }
171 }
172 }
173
174 const char* path = G4FindDataDir("G4LEDATA");
175 if (!path)
176 {
177 G4Exception("G4BremsstrahlungParameters::LoadData",
178 "em0006",FatalException,"G4LEDATA environment variable not set");
179 return;
180 }
181
182 G4String pathString(path);
183 G4String path2("/ioni/ion-sp-");
184 pathString += path2;
185
186 G4double energy, sum;
187
188 std::size_t nZ = activeZ.size();
189
190 for (std::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");
203 G4Exception("G4eIonisationParameters::LoadData",
204 "em0003",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<G4VEMDataSet*> p;
219 for (std::size_t k=0; k<length; ++k)
220 {
221 G4VDataSetAlgorithm* inter = new G4LinLogLogInterpolation();
222 G4VEMDataSet* composite = new G4CompositeEMDataSet(inter,1.,1.);
223 p.push_back(composite);
224 }
225
226 G4int shell = 0;
227 std::vector<G4DataVector*> a;
228 a.resize(length);
229 G4DataVector e;
230 G4bool isReady = false;
231 do {
232 file >> energy >> sum;
233 //Check if energy is valid
234 if (energy < -2)
235 {
236 G4String excep("invalid data file");
237 G4Exception("G4eIonisationParameters::LoadData",
238 "em0005",FatalException,excep);
239 return;
240 }
241
242 if (energy == -2) break;
243
244 if (energy > -1) {
245 // new energy
246 if(!isReady) {
247 isReady = true;
248 e.clear();
249 for (std::size_t j=0; j<length; ++j) {
250 a[j] = new G4DataVector();
251 }
252 }
253 e.push_back(energy);
254 a[0]->push_back(sum);
255 for (std::size_t j=1; j<length; ++j) {
256 G4double qRead;
257 file >> qRead;
258 a[j]->push_back(qRead);
259 }
260
261 } else {
262
263 // End of set for a shell, fill the map
264 for (std::size_t k=0; k<length; ++k) {
265
266 G4VDataSetAlgorithm* interp;
267 if(0 == k) { interp = new G4LinLogLogInterpolation(); }
268 else { interp = new G4LogLogInterpolation(); }
269
270 G4DataVector* eVector = new G4DataVector;
271 std::size_t eSize = e.size();
272 for (std::size_t sLocal=0; sLocal<eSize; ++sLocal) {
273 eVector->push_back(e[sLocal]);
274 }
275 G4VEMDataSet* set = new G4EMDataSet(shell,eVector,a[k],interp,1.,1.);
276
277 p[k]->AddComponent(set);
278 }
279
280 // next shell
281 ++shell;
282 isReady = false;
283 }
284 } while (energy > -2);
285
286 file.close();
287
288 for (G4int kk=0; kk<(G4int)length; ++kk)
289 {
290 G4int id = Z*100 + kk;
291 param[id] = p[kk];
292 }
293 }
294
295 G4String pathString_a(path);
296 G4String name_a = pathString_a + G4String("/ioni/ion-ex-av.dat");
297 std::ifstream file_a(name_a);
298 std::filebuf* lsdp_a = file_a.rdbuf();
299 G4String pathString_b(path);
300 G4String name_b = pathString_b + G4String("/ioni/ion-ex-sig.dat");
301 std::ifstream file_b(name_b);
302 std::filebuf* lsdp_b = file_b.rdbuf();
303
304 if (! (lsdp_a->is_open()) ) {
305 G4String excep = G4String("cannot open file ")
306 + name_a;
307 G4Exception("G4eIonisationParameters::LoadData",
308 "em0003",FatalException,excep);
309 }
310 if (! (lsdp_b->is_open()) ) {
311 G4String excep = G4String("cannot open file ")
312 + name_b;
313 G4Exception("G4eIonisationParameters::LoadData",
314 "em0003",FatalException,excep);
315 }
316
317 // The file is organized into two columns:
318 // 1st column is the energy
319 // 2nd column is the corresponding value
320 // The file terminates with the pattern: -1 -1
321 // -2 -2
322
323 G4double ener, ener1, sig, sig1;
324 G4int z = 0;
325
326 G4DataVector e;
327 e.clear();
328 G4DataVector d;
329 d.clear();
330
331 do {
332 file_a >> ener >> sig;
333 file_b >> ener1 >> sig1;
334
335 if(ener != ener1) {
336 G4cout << "G4eIonisationParameters: problem in excitation data "
337 << "ener= " << ener
338 << " ener1= " << ener1
339 << G4endl;
340 }
341
342 // End of file
343 if (ener == -2) {
344 break;
345
346 // End of next element
347 } else if (ener == -1) {
348
349 z++;
350 G4double Z = (G4double)z;
351
352 // fill map if Z is used
353 if (activeZ.contains(Z)) {
354
355 G4VDataSetAlgorithm* inter = new G4LinInterpolation();
356 G4DataVector* eVector = new G4DataVector;
357 G4DataVector* dVector = new G4DataVector;
358 std::size_t eSize = e.size();
359 for (std::size_t sLocal2=0; sLocal2<eSize; ++sLocal2) {
360 eVector->push_back(e[sLocal2]);
361 dVector->push_back(d[sLocal2]);
362 }
363 G4VEMDataSet* set = new G4EMDataSet(z,eVector,dVector,inter,1.,1.);
364 excit[z] = set;
365 }
366 e.clear();
367 d.clear();
368
369 } else {
370
371 e.push_back(ener*MeV);
372 d.push_back(sig1*sig*barn*MeV);
373 }
374 } while (ener != -2);
375
376 file_a.close();
377
378 }
379
380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
381
382 void G4eIonisationParameters::PrintData() const
383 {
384 G4cout << G4endl;
385 G4cout << "===== G4eIonisationParameters =====" << G4endl;
386 G4cout << G4endl;
387
388 std::size_t nZ = activeZ.size();
389 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
390
391 for (std::size_t i=0; i<nZ; i++) {
392 G4int Z = (G4int)activeZ[i];
393
394 for (G4int j=0; j<(G4int)length; ++j) {
395
396 G4int index = Z*100 + j;
397
398 pos = param.find(index);
399 if (pos!= param.cend()) {
400 G4VEMDataSet* dataSet = (*pos).second;
401 std::size_t nShells = dataSet->NumberOfComponents();
402
403 for (G4int k=0; k<(G4int)nShells; ++k) {
404
405 G4cout << "===== Z= " << Z << " shell= " << k
406 << " parameter[" << j << "] ====="
407 << G4endl;
408 const G4VEMDataSet* comp = dataSet->GetComponent(k);
409 comp->PrintData();
410 }
411 }
412 }
413 }
414 G4cout << "====================================" << G4endl;
415 }
416
417
418