| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //
28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@
29 //
30 // History:
31 // -----------
32 // 31 Jul 2001 MGP Created, with dumm
33 // 12.09.01 V.Ivanchenko Add param and inte
34 // 04.10.01 V.Ivanchenko Add BindingEnergy
35 // 25.10.01 MGP Many bug fixes, mo
36 // management of poin
37 // 29.11.01 V.Ivanchenko New parametrisatio
38 // 30.05.02 V.Ivanchenko Format and names o
39 // chenged to "ion-..
40 // 17.02.04 V.Ivanchenko Increase buffer si
41 //
42 // -------------------------------------------
43
44 #include <fstream>
45 #include <sstream>
46
47 #include "G4RDeIonisationParameters.hh"
48 #include "G4RDVEMDataSet.hh"
49 #include "G4RDShellEMDataSet.hh"
50 #include "G4RDEMDataSet.hh"
51 #include "G4RDCompositeEMDataSet.hh"
52 #include "G4RDLinInterpolation.hh"
53 #include "G4RDLogLogInterpolation.hh"
54 #include "G4RDLinLogLogInterpolation.hh"
55 #include "G4RDSemiLogInterpolation.hh"
56 #include "G4SystemOfUnits.hh"
57 #include "G4Material.hh"
58 #include "G4DataVector.hh"
59
60 G4RDeIonisationParameters:: G4RDeIonisationPar
61 : zMin(minZ), zMax(maxZ),
62 length(24)
63 {
64 LoadData();
65 }
66
67
68 G4RDeIonisationParameters::~G4RDeIonisationPar
69 {
70 // Reset the map of data sets: remove the da
71 std::map<G4int,G4RDVEMDataSet*,std::less<G4i
72
73 for (pos = param.begin(); pos != param.end()
74 {
75 G4RDVEMDataSet* dataSet = (*pos).second;
76 delete dataSet;
77 }
78
79 for (pos = excit.begin(); pos != excit.end()
80 {
81 G4RDVEMDataSet* dataSet = (*pos).second;
82 delete dataSet;
83 }
84
85 activeZ.clear();
86 }
87
88
89 G4double G4RDeIonisationParameters::Parameter(
90
91
92 {
93 G4double value = 0.;
94 G4int id = Z*100 + parameterIndex;
95 std::map<G4int,G4RDVEMDataSet*,std::less<G4i
96
97 pos = param.find(id);
98 if (pos!= param.end()) {
99 G4RDVEMDataSet* dataSet = (*pos).second;
100 G4int nShells = dataSet->NumberOfComponent
101
102 if(shellIndex < nShells) {
103 const G4RDVEMDataSet* component = dataSe
104 const G4DataVector ener = component->Get
105 G4double ee = std::max(ener.front(),std:
106 value = component->FindValue(ee);
107 } else {
108 G4cout << "WARNING: G4IonisationParamete
109 << "has no parameters for shell=
110 << "; Z= " << Z
111 << G4endl;
112 }
113 } else {
114 G4cout << "WARNING: G4IonisationParameters
115 << "did not find ID = "
116 << shellIndex << G4endl;
117 }
118
119 return value;
120 }
121
122 G4double G4RDeIonisationParameters::Excitation
123 {
124 G4double value = 0.;
125 std::map<G4int,G4RDVEMDataSet*,std::less<G4i
126
127 pos = excit.find(Z);
128 if (pos!= excit.end()) {
129 G4RDVEMDataSet* dataSet = (*pos).second;
130
131 const G4DataVector ener = dataSet->GetEner
132 G4double ee = std::max(ener.front(),std::m
133 value = dataSet->FindValue(ee);
134 } else {
135 G4cout << "WARNING: G4IonisationParameters
136 << "did not find ID = "
137 << Z << G4endl;
138 }
139
140 return value;
141 }
142
143
144 void G4RDeIonisationParameters::LoadData()
145 {
146 // ---------------------------------------
147 // Please document what are the parameters
148 // ---------------------------------------
149
150 // define active elements
151
152 const G4MaterialTable* materialTable = G4Mat
153 if (materialTable == 0)
154 G4Exception("G4RDeIonisationParameters::L
155 FatalException, "No MaterialT
156
157 G4int nMaterials = G4Material::GetNumberOfMa
158
159 for (G4int m=0; m<nMaterials; m++) {
160
161 const G4Material* material= (*materialTabl
162 const G4ElementVector* elementVector = mat
163 const size_t nElements = material->GetNumb
164
165 for (size_t iEl=0; iEl<nElements; iEl++) {
166 G4Element* element = (*elementVector)[iE
167 G4double Z = element->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 G4String excep("G4LEDATA environment var
178 G4Exception("G4RDeIonisationParameters::
179 FatalException, 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("Data file: ")
202 + name + G4String(" not found. The versi
203 G4Exception("G4RDeIonisationParameters::
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 t
212 // 2d line:
213 // 3 energy (MeV) subdividing different ap
214 // 20 point on the spectrum
215 // The shell terminates with the pattern:
216 // The file terminates with the pattern: -
217
218 std::vector<G4RDVEMDataSet*> p;
219 for (size_t k=0; k<length; k++)
220 {
221 G4RDVDataSetAlgorithm* inter = new G4RDLinLo
222 G4RDVEMDataSet* composite = new G4RDComposit
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 ma
251 for (size_t k=0; k<length; k++) {
252
253 G4RDVDataSetAlgorithm* interp;
254 if(0 == k) interp = new G4RDLinLogL
255 else interp = new G4RDLogLogI
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(sh
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("/
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("/
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 fi
296 + name_a;
297 G4Exception("G4RDeIonisationParameters::L
298 FatalException, excep);
299 }
300 if (! (lsdp_b->is_open()) ) {
301 G4String excep = G4String("Cannot open fi
302 + name_b;
303 G4Exception("G4RDeIonisationParameters::L
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
311 // -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: pr
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
345 G4RDVDataSetAlgorithm* inter = new G4RDLinI
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,eV
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
371 void G4RDeIonisationParameters::PrintData() co
372 {
373 G4cout << G4endl;
374 G4cout << "===== G4RDeIonisationParameters =
375 G4cout << G4endl;
376
377 size_t nZ = activeZ.size();
378 std::map<G4int,G4RDVEMDataSet*,std::less<G4i
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).secon
390 size_t nShells = dataSet->NumberOfComp
391
392 for (size_t k=0; k<nShells; k++) {
393
394 G4cout << "===== Z= " << Z << " shel
395 << " parameter[" << j << "]
396 << G4endl;
397 const G4RDVEMDataSet* comp = dataSet
398 comp->PrintData();
399 }
400 }
401 }
402 }
403 G4cout << "=================================
404 }
405
406
407