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

Last change on this file since 1197 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
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.25 2009/06/10 13:32:36 mantero Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
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// 23.03.09 L.Pandola Updated warning message
44//
45// -------------------------------------------------------------------
46
47#include "G4eIonisationParameters.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
62G4eIonisationParameters:: G4eIonisationParameters(G4int minZ, G4int maxZ)
63 : zMin(minZ), zMax(maxZ),
64 length(24)
65{
66 LoadData();
67}
68
69
70G4eIonisationParameters::~G4eIonisationParameters()
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
91G4double G4eIonisationParameters::Parameter(G4int Z, G4int shellIndex,
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
124G4double G4eIonisationParameters::Excitation(G4int Z, G4double e) const
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
146void G4eIonisationParameters::LoadData()
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: no MaterialTable found)");
157
158 G4int nMaterials = G4Material::GetNumberOfMaterials();
159
160 for (G4int m=0; m<nMaterials; m++) {
161
162 const G4Material* material= (*materialTable)[m];
163 const G4ElementVector* elementVector = material->GetElementVector();
164 const size_t nElements = material->GetNumberOfElements();
165
166 for (size_t iEl=0; iEl<nElements; iEl++) {
167 G4Element* element = (*elementVector)[iEl];
168 G4double Z = element->GetZ();
169 if (!(activeZ.contains(Z))) {
170 activeZ.push_back(Z);
171 }
172 }
173 }
174
175 char* path = getenv("G4LEDATA");
176 if (!path)
177 {
178 G4String excep("G4eIonisationParameters - G4LEDATA environment variable not set");
179 G4Exception(excep);
180 }
181
182 G4String pathString(path);
183 G4String path2("/ioni/ion-sp-");
184 pathString += path2;
185
186 G4double energy, sum;
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("G4IonisationParameters - data file: ")
202 + name + G4String(" not found. Please check and/or update G4LEDATA");
203 G4Exception(excep);
204 }
205
206 // The file is organized into:
207 // For each shell there are two lines:
208 // 1st line:
209 // 1st column is the energy of incident e-,
210 // 2d column is the parameter of screan term;
211 // 2d line:
212 // 3 energy (MeV) subdividing different approximation area of the spectrum
213 // 20 point on the spectrum
214 // The shell terminates with the pattern: -1 -1
215 // The file terminates with the pattern: -2 -2
216
217 std::vector<G4VEMDataSet*> p;
218 for (size_t k=0; k<length; k++)
219 {
220 G4VDataSetAlgorithm* inter = new G4LinLogLogInterpolation();
221 G4VEMDataSet* composite = new G4CompositeEMDataSet(inter,1.,1.);
222 p.push_back(composite);
223 }
224
225 G4int shell = 0;
226 std::vector<G4DataVector*> a;
227 for (size_t j=0; j<length; j++)
228 {
229 G4DataVector* aa = new G4DataVector();
230 a.push_back(aa);
231 }
232 G4DataVector e;
233 e.clear();
234 do {
235 file >> energy >> sum;
236 if (energy == -2) break;
237
238 if (energy > -1) {
239 e.push_back(energy);
240 a[0]->push_back(sum);
241 for (size_t j=0; j<length-1; j++) {
242 G4double qRead;
243 file >> qRead;
244 a[j + 1]->push_back(qRead);
245 }
246
247 } else {
248
249 // End of set for a shell, fill the map
250 for (size_t k=0; k<length; k++) {
251
252 G4VDataSetAlgorithm* interp;
253 if(0 == k) interp = new G4LinLogLogInterpolation();
254 else interp = new G4LogLogInterpolation();
255
256 G4DataVector* eVector = new G4DataVector;
257 size_t eSize = e.size();
258 for (size_t s=0; s<eSize; s++) {
259 eVector->push_back(e[s]);
260 }
261 G4VEMDataSet* set = new G4EMDataSet(shell,eVector,a[k],interp,1.,1.);
262
263 p[k]->AddComponent(set);
264 }
265
266 // clear vectors
267 for (size_t j2=0; j2<length; j2++) {
268 a[j2] = new G4DataVector();
269 }
270 shell++;
271 e.clear();
272 }
273 } while (energy > -2);
274
275 file.close();
276
277 for (size_t kk=0; kk<length; kk++)
278 {
279 G4int id = Z*100 + kk;
280 param[id] = p[kk];
281 }
282 }
283
284 G4String pathString_a(path);
285 G4String name_a = pathString_a + G4String("/ioni/ion-ex-av.dat");
286 std::ifstream file_a(name_a);
287 std::filebuf* lsdp_a = file_a.rdbuf();
288 G4String pathString_b(path);
289 G4String name_b = pathString_b + G4String("/ioni/ion-ex-sig.dat");
290 std::ifstream file_b(name_b);
291 std::filebuf* lsdp_b = file_b.rdbuf();
292
293 if (! (lsdp_a->is_open()) ) {
294 G4String excep = G4String("G4eIonisationParameters: cannot open file ")
295 + name_a;
296 G4Exception(excep);
297 }
298 if (! (lsdp_b->is_open()) ) {
299 G4String excep = G4String("G4eIonisationParameters: cannot open file ")
300 + name_b;
301 G4Exception(excep);
302 }
303
304 // The file is organized into two columns:
305 // 1st column is the energy
306 // 2nd column is the corresponding value
307 // The file terminates with the pattern: -1 -1
308 // -2 -2
309
310 G4double ener, ener1, sig, sig1;
311 G4int z = 0;
312
313 G4DataVector e;
314 e.clear();
315 G4DataVector d;
316 d.clear();
317
318 do {
319 file_a >> ener >> sig;
320 file_b >> ener1 >> sig1;
321
322 if(ener != ener1) {
323 G4cout << "G4eIonisationParameters: problem in excitation data "
324 << "ener= " << ener
325 << " ener1= " << ener1
326 << G4endl;
327 }
328
329 // End of file
330 if (ener == -2) {
331 break;
332
333 // End of next element
334 } else if (ener == -1) {
335
336 z++;
337 G4double Z = (G4double)z;
338
339 // fill map if Z is used
340 if (activeZ.contains(Z)) {
341
342 G4VDataSetAlgorithm* inter = new G4LinInterpolation();
343 G4DataVector* eVector = new G4DataVector;
344 G4DataVector* dVector = new G4DataVector;
345 size_t eSize = e.size();
346 for (size_t s=0; s<eSize; s++) {
347 eVector->push_back(e[s]);
348 dVector->push_back(d[s]);
349 }
350 G4VEMDataSet* set = new G4EMDataSet(z,eVector,dVector,inter,1.,1.);
351 excit[z] = set;
352 }
353 e.clear();
354 d.clear();
355
356 } else {
357
358 e.push_back(ener*MeV);
359 d.push_back(sig1*sig*barn*MeV);
360 }
361 } while (ener != -2);
362
363 file_a.close();
364
365}
366
367
368void G4eIonisationParameters::PrintData() const
369{
370 G4cout << G4endl;
371 G4cout << "===== G4eIonisationParameters =====" << G4endl;
372 G4cout << G4endl;
373
374 size_t nZ = activeZ.size();
375 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
376
377 for (size_t i=0; i<nZ; i++) {
378 G4int Z = (G4int)activeZ[i];
379
380 for (size_t j=0; j<length; j++) {
381
382 G4int index = Z*100 + j;
383
384 pos = param.find(index);
385 if (pos!= param.end()) {
386 G4VEMDataSet* dataSet = (*pos).second;
387 size_t nShells = dataSet->NumberOfComponents();
388
389 for (size_t k=0; k<nShells; k++) {
390
391 G4cout << "===== Z= " << Z << " shell= " << k
392 << " parameter[" << j << "] ====="
393 << G4endl;
394 const G4VEMDataSet* comp = dataSet->GetComponent(k);
395 comp->PrintData();
396 }
397 }
398 }
399 }
400 G4cout << "====================================" << G4endl;
401}
402
403
Note: See TracBrowser for help on using the repository browser.