source: trunk/source/processes/electromagnetic/lowenergy/src/G4eIonisationParameters.cc@ 1005

Last change on this file since 1005 was 991, checked in by garnier, 17 years ago

update
File size: 11.6 KB
Line 
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,v 1.23 2006/06/29 19:42:02 gunter Exp $
28// GEANT4 tag $Name: geant4-09-02 $
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 "G4eIonisationParameters.hh"
47#include "G4VEMDataSet.hh"
48#include "G4ShellEMDataSet.hh"
49#include "G4EMDataSet.hh"
50#include "G4CompositeEMDataSet.hh"
51#include "G4LinInterpolation.hh"
52#include "G4LogLogInterpolation.hh"
53#include "G4LinLogLogInterpolation.hh"
54#include "G4SemiLogInterpolation.hh"
55#include "G4Material.hh"
56#include "G4DataVector.hh"
57#include <fstream>
58#include <sstream>
59
60
61G4eIonisationParameters:: G4eIonisationParameters(G4int minZ, G4int maxZ)
62 : zMin(minZ), zMax(maxZ),
63 length(24)
64{
65 LoadData();
66}
67
68
69G4eIonisationParameters::~G4eIonisationParameters()
70{
71 // Reset the map of data sets: remove the data sets from the map
72 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::iterator pos;
73
74 for (pos = param.begin(); pos != param.end(); ++pos)
75 {
76 G4VEMDataSet* dataSet = (*pos).second;
77 delete dataSet;
78 }
79
80 for (pos = excit.begin(); pos != excit.end(); ++pos)
81 {
82 G4VEMDataSet* dataSet = (*pos).second;
83 delete dataSet;
84 }
85
86 activeZ.clear();
87}
88
89
90G4double 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 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
97
98 pos = param.find(id);
99 if (pos!= param.end()) {
100 G4VEMDataSet* dataSet = (*pos).second;
101 G4int nShells = dataSet->NumberOfComponents();
102
103 if(shellIndex < nShells) {
104 const G4VEMDataSet* component = dataSet->GetComponent(shellIndex);
105 const G4DataVector ener = component->GetEnergies(0);
106 G4double ee = std::max(ener.front(),std::min(ener.back(),e));
107 value = component->FindValue(ee);
108 } else {
109 G4cout << "WARNING: G4IonisationParameters::FindParameter "
110 << "has no parameters for shell= " << shellIndex
111 << "; Z= " << Z
112 << G4endl;
113 }
114 } else {
115 G4cout << "WARNING: G4IonisationParameters::Parameter "
116 << "did not find ID = "
117 << shellIndex << G4endl;
118 }
119
120 return value;
121}
122
123G4double G4eIonisationParameters::Excitation(G4int Z, G4double e) const
124{
125 G4double value = 0.;
126 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
127
128 pos = excit.find(Z);
129 if (pos!= excit.end()) {
130 G4VEMDataSet* dataSet = (*pos).second;
131
132 const G4DataVector ener = dataSet->GetEnergies(0);
133 G4double ee = std::max(ener.front(),std::min(ener.back(),e));
134 value = dataSet->FindValue(ee);
135 } else {
136 G4cout << "WARNING: G4IonisationParameters::Excitation "
137 << "did not find ID = "
138 << Z << G4endl;
139 }
140
141 return value;
142}
143
144
145void 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 == 0)
155 G4Exception("G4eIonisationParameters: no MaterialTable found)");
156
157 G4int nMaterials = G4Material::GetNumberOfMaterials();
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 = getenv("G4LEDATA");
175 if (!path)
176 {
177 G4String excep("G4eIonisationParameters - G4LEDATA environment variable not set");
178 G4Exception(excep);
179 }
180
181 G4String pathString(path);
182 G4String path2("/ioni/ion-sp-");
183 pathString += path2;
184
185 G4double energy, sum;
186
187 size_t nZ = activeZ.size();
188
189 for (size_t i=0; i<nZ; i++) {
190
191 G4int Z = (G4int)activeZ[i];
192 std::ostringstream ost;
193 ost << pathString << Z << ".dat";
194 G4String name(ost.str());
195
196 std::ifstream file(name);
197 std::filebuf* lsdp = file.rdbuf();
198
199 if (! (lsdp->is_open()) ) {
200 G4String excep = G4String("G4IonisationParameters - data file: ")
201 + name + G4String(" not found. The version 1.# of G4LEDATA should be used");
202 G4Exception(excep);
203 }
204
205 // The file is organized into:
206 // For each shell there are two lines:
207 // 1st line:
208 // 1st column is the energy of incident e-,
209 // 2d column is the parameter of screan term;
210 // 2d line:
211 // 3 energy (MeV) subdividing different approximation area of the spectrum
212 // 20 point on the spectrum
213 // The shell terminates with the pattern: -1 -1
214 // The file terminates with the pattern: -2 -2
215
216 std::vector<G4VEMDataSet*> p;
217 for (size_t k=0; k<length; k++)
218 {
219 G4VDataSetAlgorithm* inter = new G4LinLogLogInterpolation();
220 G4VEMDataSet* composite = new G4CompositeEMDataSet(inter,1.,1.);
221 p.push_back(composite);
222 }
223
224 G4int shell = 0;
225 std::vector<G4DataVector*> a;
226 for (size_t j=0; j<length; j++)
227 {
228 G4DataVector* aa = new G4DataVector();
229 a.push_back(aa);
230 }
231 G4DataVector e;
232 e.clear();
233 do {
234 file >> energy >> sum;
235 if (energy == -2) break;
236
237 if (energy > -1) {
238 e.push_back(energy);
239 a[0]->push_back(sum);
240 for (size_t j=0; j<length-1; j++) {
241 G4double qRead;
242 file >> qRead;
243 a[j + 1]->push_back(qRead);
244 }
245
246 } else {
247
248 // End of set for a shell, fill the map
249 for (size_t k=0; k<length; k++) {
250
251 G4VDataSetAlgorithm* interp;
252 if(0 == k) interp = new G4LinLogLogInterpolation();
253 else interp = new G4LogLogInterpolation();
254
255 G4DataVector* eVector = new G4DataVector;
256 size_t eSize = e.size();
257 for (size_t s=0; s<eSize; s++) {
258 eVector->push_back(e[s]);
259 }
260 G4VEMDataSet* set = new G4EMDataSet(shell,eVector,a[k],interp,1.,1.);
261
262 p[k]->AddComponent(set);
263 }
264
265 // clear vectors
266 for (size_t j2=0; j2<length; j2++) {
267 a[j2] = new G4DataVector();
268 }
269 shell++;
270 e.clear();
271 }
272 } while (energy > -2);
273
274 file.close();
275
276 for (size_t kk=0; kk<length; kk++)
277 {
278 G4int id = Z*100 + kk;
279 param[id] = p[kk];
280 }
281 }
282
283 G4String pathString_a(path);
284 G4String name_a = pathString_a + G4String("/ioni/ion-ex-av.dat");
285 std::ifstream file_a(name_a);
286 std::filebuf* lsdp_a = file_a.rdbuf();
287 G4String pathString_b(path);
288 G4String name_b = pathString_b + G4String("/ioni/ion-ex-sig.dat");
289 std::ifstream file_b(name_b);
290 std::filebuf* lsdp_b = file_b.rdbuf();
291
292 if (! (lsdp_a->is_open()) ) {
293 G4String excep = G4String("G4eIonisationParameters: cannot open file ")
294 + name_a;
295 G4Exception(excep);
296 }
297 if (! (lsdp_b->is_open()) ) {
298 G4String excep = G4String("G4eIonisationParameters: cannot open file ")
299 + name_b;
300 G4Exception(excep);
301 }
302
303 // The file is organized into two columns:
304 // 1st column is the energy
305 // 2nd column is the corresponding value
306 // The file terminates with the pattern: -1 -1
307 // -2 -2
308
309 G4double ener, ener1, sig, sig1;
310 G4int z = 0;
311
312 G4DataVector e;
313 e.clear();
314 G4DataVector d;
315 d.clear();
316
317 do {
318 file_a >> ener >> sig;
319 file_b >> ener1 >> sig1;
320
321 if(ener != ener1) {
322 G4cout << "G4eIonisationParameters: problem in excitation data "
323 << "ener= " << ener
324 << " ener1= " << ener1
325 << G4endl;
326 }
327
328 // End of file
329 if (ener == -2) {
330 break;
331
332 // End of next element
333 } else if (ener == -1) {
334
335 z++;
336 G4double Z = (G4double)z;
337
338 // fill map if Z is used
339 if (activeZ.contains(Z)) {
340
341 G4VDataSetAlgorithm* inter = new G4LinInterpolation();
342 G4DataVector* eVector = new G4DataVector;
343 G4DataVector* dVector = new G4DataVector;
344 size_t eSize = e.size();
345 for (size_t s=0; s<eSize; s++) {
346 eVector->push_back(e[s]);
347 dVector->push_back(d[s]);
348 }
349 G4VEMDataSet* set = new G4EMDataSet(z,eVector,dVector,inter,1.,1.);
350 excit[z] = set;
351 }
352 e.clear();
353 d.clear();
354
355 } else {
356
357 e.push_back(ener*MeV);
358 d.push_back(sig1*sig*barn*MeV);
359 }
360 } while (ener != -2);
361
362 file_a.close();
363
364}
365
366
367void G4eIonisationParameters::PrintData() const
368{
369 G4cout << G4endl;
370 G4cout << "===== G4eIonisationParameters =====" << G4endl;
371 G4cout << G4endl;
372
373 size_t nZ = activeZ.size();
374 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
375
376 for (size_t i=0; i<nZ; i++) {
377 G4int Z = (G4int)activeZ[i];
378
379 for (size_t j=0; j<length; j++) {
380
381 G4int index = Z*100 + j;
382
383 pos = param.find(index);
384 if (pos!= param.end()) {
385 G4VEMDataSet* dataSet = (*pos).second;
386 size_t nShells = dataSet->NumberOfComponents();
387
388 for (size_t k=0; k<nShells; k++) {
389
390 G4cout << "===== Z= " << Z << " shell= " << k
391 << " parameter[" << j << "] ====="
392 << G4endl;
393 const G4VEMDataSet* comp = dataSet->GetComponent(k);
394 comp->PrintData();
395 }
396 }
397 }
398 }
399 G4cout << "====================================" << G4endl;
400}
401
402
Note: See TracBrowser for help on using the repository browser.