| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // 23 //
>> 24 // $Id: G4ShellEMDataSet.cc,v 1.10 2003/06/16 17:00:26 gunter Exp $
>> 25 // GEANT4 tag $Name: geant4-06-00 $
27 // 26 //
28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@ 27 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
29 // 28 //
30 // History: 29 // History:
31 // ----------- 30 // -----------
32 // 1 Aug 2001 MGP Created << 31 // 1 Aug 2001 MGP Created
33 // << 32 // 09.10.01 V.Ivanchenko Add case z=0
34 // 09.10.01 V.Ivanchenko Add case z=0 <<
35 // <<
36 // 9 Mar 2008 MGP Cleaned up unread <<
37 // (Further clean-up <<
38 // <<
39 // 15 Jul 2009 Nicolas A. Karakatsanis <<
40 // <<
41 // - LoadNonLogData <<
42 // dataset. It is <<
43 // <<
44 // - LoadData method <<
45 // It retrieves th <<
46 // respective log <<
47 // <<
48 // - SetLogEnergiesD <<
49 // The EM data set <<
50 // These initializ <<
51 // operations <<
52 // <<
53 // 33 //
54 // ------------------------------------------- 34 // -------------------------------------------------------------------
55 35
56 #include "G4ShellEMDataSet.hh" 36 #include "G4ShellEMDataSet.hh"
57 #include "G4EMDataSet.hh" 37 #include "G4EMDataSet.hh"
58 #include "G4VDataSetAlgorithm.hh" 38 #include "G4VDataSetAlgorithm.hh"
59 #include <fstream> 39 #include <fstream>
60 #include <sstream> << 40 #include <strstream>
61 41
62 //....oooOO0OOooo........oooOO0OOooo........oo <<
63 G4ShellEMDataSet::G4ShellEMDataSet(G4int zeta, <<
64 G4double eUnit, <<
65 G4double dataUnit) <<
66 : <<
67 algorithm(algo), <<
68 unitEnergies(eUnit), <<
69 unitData(dataUnit), <<
70 z(zeta) <<
71 { <<
72 if (algorithm == nullptr) <<
73 G4Exception("G4ShellEMDataSet::G4ShellEMDa <<
74 FatalErrorInArgument, "Interpolation == 0" <<
75 } <<
76 42
77 //....oooOO0OOooo........oooOO0OOooo........oo << 43 G4ShellEMDataSet::G4ShellEMDataSet(G4int Z,
78 G4ShellEMDataSet::~G4ShellEMDataSet() << 44 const G4VDataSetAlgorithm* interpolation,
>> 45 G4double unitE, G4double unitData)
>> 46 :z(Z), algorithm(interpolation)
79 { 47 {
80 CleanUpComponents(); << 48 nComponents = 0;
81 if (algorithm) delete algorithm; << 49 unit1 = unitE;
>> 50 unit2 = unitData;
82 } 51 }
83 52
84 //....oooOO0OOooo........oooOO0OOooo........oo << 53 G4ShellEMDataSet::G4ShellEMDataSet(G4int Z, const G4String& dataFile,
85 G4double G4ShellEMDataSet::FindValue(G4double << 54 const G4VDataSetAlgorithm* interpolation,
>> 55 G4double unitE, G4double unitData)
>> 56 :z(Z), algorithm(interpolation)
86 { 57 {
87 // Returns the sum over the shells correspon << 58 nComponents = 0;
88 G4double value = 0.; << 59 unit1 = unitE;
89 << 60 unit2 = unitData;
90 std::vector<G4VEMDataSet *>::const_iterator << 61 LoadData(dataFile);
91 std::vector<G4VEMDataSet *>::const_iterator <<
92 <<
93 while (i != end) <<
94 { <<
95 value += (*i)->FindValue(energy); <<
96 i++; <<
97 } <<
98 <<
99 return value; <<
100 } 62 }
101 63
102 //....oooOO0OOooo........oooOO0OOooo........oo << 64 G4ShellEMDataSet::~G4ShellEMDataSet()
103 << 65 {
104 void G4ShellEMDataSet::PrintData() const << 66 for (size_t i=0; i<nComponents; i++)
105 { <<
106 const G4int n = (G4int)NumberOfComponents(); <<
107 <<
108 G4cout << "The data set has " << n << " comp <<
109 G4cout << G4endl; <<
110 <<
111 G4int i = 0; <<
112 <<
113 while (i < n) <<
114 { 67 {
115 G4cout << "--- Component " << i << " --- << 68 G4VEMDataSet* dataSet = components[i];
116 GetComponent(i)->PrintData(); << 69 delete dataSet;
117 ++i; <<
118 } 70 }
>> 71 delete algorithm;
119 } 72 }
120 73
121 //....oooOO0OOooo........oooOO0OOooo........oo << 74 G4double G4ShellEMDataSet::FindValue(G4double e, G4int ) const
122 <<
123 void G4ShellEMDataSet::SetEnergiesData(G4DataV <<
124 G4DataVector* data, <<
125 G4int componentId) <<
126 { 75 {
127 G4VEMDataSet* component = components[compone << 76 // Returns the sum over the shells corresponding to e
128 if (component) << 77 G4double value = 0.;
129 { <<
130 component->SetEnergiesData(energies, dat <<
131 return; <<
132 } <<
133 <<
134 G4String msg = "component "; <<
135 msg += static_cast<char>(componentId); <<
136 msg += " not found"; <<
137 <<
138 G4Exception("G4ShellEMDataSet::SetEnergiesDa <<
139 } <<
140 <<
141 //....oooOO0OOooo........oooOO0OOooo........oo <<
142 78
143 void G4ShellEMDataSet::SetLogEnergiesData(G4Da << 79 for (size_t i=0; i<nComponents; i++)
144 G4Da <<
145 G4Da <<
146 G4Da <<
147 G4in <<
148 { <<
149 G4VEMDataSet* component = components[compone <<
150 if (component) <<
151 { 80 {
152 component->SetLogEnergiesData(energies, << 81 G4VEMDataSet* component = components[i];
153 return; << 82 G4double shellValue = component->FindValue(e);
>> 83 value = value + shellValue;
154 } 84 }
155 85
156 G4String msg = "component "; << 86 return value;
157 msg += static_cast<char>(componentId); <<
158 msg += " not found"; <<
159 G4Exception("G4ShellEMDataSet::SetLogEnergie <<
160 FatalErrorInArgument ,msg); <<
161 } 87 }
162 88
163 //....oooOO0OOooo........oooOO0OOooo........oo << 89 void G4ShellEMDataSet::PrintData() const
164 <<
165 G4bool G4ShellEMDataSet::LoadData(const G4Stri <<
166 { 90 {
167 CleanUpComponents(); << 91 G4cout << "The data set has " << nComponents << " components" << G4endl;
168 92
169 G4String fullFileName = FullFileName(file); << 93 for (size_t i=0; i<nComponents; i++)
170 std::ifstream in(fullFileName); << 94 {
>> 95 G4cout << "--- Component " << i << " ---" << G4endl;
>> 96 G4VEMDataSet* component = components[i];
>> 97 component->PrintData();
>> 98 }
>> 99 }
171 100
172 if (!in.is_open()) << 101 void G4ShellEMDataSet::LoadData(const G4String& fileName)
173 { << 102 {
174 G4String message("Data file \""); << 103 // Build the complete string identifying the file with the data set
175 message += fullFileName; <<
176 message += "\" not found"; <<
177 G4Exception("G4ShellEMDataSet::LoadData( <<
178 return 0; <<
179 } <<
180 104
181 G4DataVector* orig_shell_energies = 0; << 105 char nameChar[100] = {""};
182 G4DataVector* orig_shell_data = 0; << 106 std::ostrstream ost(nameChar, 100, std::ios::out);
183 G4DataVector* log_shell_energies = 0; << 107
184 G4DataVector* log_shell_data = 0; << 108 if (z != 0) ost << fileName << z << ".dat";
185 << 109 else ost << fileName << ".dat";
186 G4double a = 0.; <<
187 G4int shellIndex = 0; <<
188 G4int k = 0; <<
189 G4int nColumns = 2; <<
190 <<
191 do <<
192 { <<
193 in >> a; <<
194 <<
195 if (a==0.) a=1e-300; <<
196 110
197 // The file is organized into four colum << 111 G4String name(nameChar);
198 // 1st column contains the values of ene << 112
199 // 2nd column contains the corresponding << 113 char* path = getenv("G4LEDATA");
200 // The file terminates with the pattern: << 114 if (!path)
201 // << 115 {
202 // << 116 G4String excep("G4ShellEMDataSet - G4LEDATA environment variable not set");
203 if (a == -1) << 117 G4Exception(excep);
204 { <<
205 if ((k%nColumns == 0) && (orig_shell_energ <<
206 { <<
207 AddComponent(new G4EMDataSet(shellIndex <<
208 orig_shell_data, log_shell_energie <<
209 log_shell_data, algorithm->Clone() <<
210 unitEnergies, unitData)); <<
211 orig_shell_energies = 0; <<
212 orig_shell_data = 0; <<
213 log_shell_energies = 0; <<
214 log_shell_data = 0; <<
215 } <<
216 } <<
217 else if (a != -2) <<
218 { <<
219 if (orig_shell_energies == 0) <<
220 { <<
221 orig_shell_energies = new G4DataVector; <<
222 orig_shell_data = new G4DataVector; <<
223 log_shell_energies = new G4DataVe <<
224 log_shell_data = new G4DataVector <<
225 } <<
226 if (k%nColumns == 0) <<
227 { <<
228 orig_shell_energies->push_back(a*unitEn <<
229 log_shell_energies->push_back(std::log1 <<
230 } <<
231 else if (k%nColumns == 1) <<
232 { <<
233 orig_shell_data->push_back(a*unitData); <<
234 log_shell_data->push_back(std::lo <<
235 } <<
236 k++; <<
237 } <<
238 else k = 1; <<
239 } 118 }
240 while (a != -2); // End of file << 119
241 << 120 G4String pathString(path);
242 delete orig_shell_energies; << 121 G4String separator("/" );
243 delete orig_shell_data; << 122 G4String dirFile = pathString + separator + name;
244 delete log_shell_energies; << 123 std::ifstream file(dirFile);
245 delete log_shell_data; << 124 std::filebuf* lsdp = file.rdbuf();
246 <<
247 return true; <<
248 } <<
249 <<
250 //....oooOO0OOooo........oooOO0OOooo........oo <<
251 <<
252 G4bool G4ShellEMDataSet::LoadNonLogData(const <<
253 { <<
254 CleanUpComponents(); <<
255 <<
256 G4String fullFileName = FullFileName(file); <<
257 std::ifstream in(fullFileName); <<
258 125
259 if (!in.is_open()) << 126 if (! (lsdp->is_open()) )
260 { 127 {
261 G4String message("G4ShellEMDataSet::Load << 128 G4String s1("G4ShellEMDataSet - data file: ");
262 message += fullFileName; << 129 G4String s2(" not found");
263 message += "\" not found"; << 130 G4String excep = s1 + dirFile + s2;
264 G4Exception("G4ShellEMDataSet::LoadNonLo << 131 G4Exception(excep);
265 return 0; <<
266 } 132 }
267 133
268 G4DataVector* orig_shell_energies = 0; << 134 G4double a = 0;
269 G4DataVector* orig_shell_data = 0; << 135 G4int k = 1;
270 << 136 G4int s = 0;
271 G4double a = 0.; << 137
272 G4int shellIndex = 0; 138 G4int shellIndex = 0;
273 G4int k = 0; << 139 G4DataVector* energies = new G4DataVector;
274 G4int nColumns = 2; << 140 G4DataVector* data = new G4DataVector;
275 141
276 do << 142 do {
277 { << 143 file >> a;
278 in >> a; << 144 G4int nColumns = 2;
279 // The file is organized into four colum << 145 if (a == -1)
280 // 1st column contains the values of ene << 146 {
281 // 2nd column contains the corresponding << 147 if (s == 0)
282 // The file terminates with the pattern: << 148 {
283 // << 149 // End of a shell data set
284 // << 150 G4VDataSetAlgorithm* algo = algorithm->Clone();
285 if (a == -1) << 151 G4VEMDataSet* dataSet = new G4EMDataSet(shellIndex,energies,data,algo);
286 { << 152 AddComponent(dataSet);
287 if ((k%nColumns == 0) && (orig_shell_energ << 153 // Start of new shell data set
288 { << 154 energies = new G4DataVector;
289 AddComponent(new G4EMDataSet(shellIndex << 155 data = new G4DataVector;
290 orig_shell_data, algorithm->Clone( << 156 shellIndex++;
291 unitEnergies, unitData)); << 157 }
292 orig_shell_energies = 0; << 158 s++;
293 orig_shell_data = 0; << 159 if (s == nColumns)
294 } <<
295 } <<
296 else if (a != -2) <<
297 { 160 {
298 if (orig_shell_energies == 0) << 161 s = 0;
299 { <<
300 orig_shell_energies = new G4DataVector; <<
301 orig_shell_data = new G4DataVector; <<
302 } <<
303 if (k%nColumns == 0) <<
304 { <<
305 orig_shell_energies->push_back(a*unitEn <<
306 } <<
307 else if (k%nColumns == 1) <<
308 { <<
309 orig_shell_data->push_back(a*unitData); <<
310 } <<
311 k++; <<
312 } 162 }
313 else k = 1; << 163 }
314 } << 164 else if (a == -2)
315 while (a != -2); // End of file << 165 {
316 << 166 // End of file; delete the empty vectors created when encountering the last -1 -1 row
317 << 167 delete energies;
318 delete orig_shell_energies; << 168 delete data;
319 delete orig_shell_data; << 169 }
320 << 170 else
321 return true; << 171 {
322 } << 172 // 1st column is energy
323 << 173 if(k%nColumns != 0)
324 //....oooOO0OOooo........oooOO0OOooo........oo << 174 {
325 << 175 G4double e = a * unit1;
326 G4bool G4ShellEMDataSet::SaveData(const G4Stri << 176 energies->push_back(e);
>> 177 k++;
>> 178 }
>> 179 else if (k%nColumns == 0)
>> 180 {
>> 181 // 2nd column is cross section
>> 182 G4double value = a * unit2;
>> 183 data->push_back(value);
>> 184 k = 1;
>> 185 }
>> 186 }
>> 187 } while (a != -2); // end of file
>> 188 file.close();
>> 189 }
>> 190
>> 191 void G4ShellEMDataSet::AddComponent(G4VEMDataSet* component)
>> 192 {
>> 193 components.push_back(component);
>> 194 nComponents++;
>> 195 }
>> 196
>> 197 const G4DataVector& G4ShellEMDataSet::GetEnergies(G4int i) const
327 { 198 {
328 G4String fullFileName = FullFileName(file); << 199 const G4VEMDataSet* component = GetComponent(i);
329 std::ofstream out(fullFileName); << 200 return (component->GetEnergies(i));
330 <<
331 if (!out.is_open()) <<
332 { <<
333 G4String message("Cannot open \""); <<
334 message += fullFileName; <<
335 message += "\""; <<
336 G4Exception("G4EMDataSet::SaveData()","e <<
337 } <<
338 <<
339 const G4int n = (G4int)NumberOfComponents(); <<
340 G4int k = 0; <<
341 <<
342 while (k < n) <<
343 { <<
344 const G4VEMDataSet* component = GetCompo <<
345 <<
346 if (component) <<
347 { <<
348 const G4DataVector& energies = component-> <<
349 const G4DataVector& data = component->GetD <<
350 <<
351 G4DataVector::const_iterator i = energies. <<
352 G4DataVector::const_iterator endI = energi <<
353 G4DataVector::const_iterator j = data.cbeg <<
354 <<
355 while (i != endI) <<
356 { <<
357 out.precision(10); <<
358 out.width(15); <<
359 out.setf(std::ofstream::left); <<
360 out << ((*i)/unitEnergies) << ' '; <<
361 <<
362 out.precision(10); <<
363 out.width(15); <<
364 out.setf(std::ofstream::left); <<
365 out << ((*j)/unitData) << std::endl; <<
366 i++; <<
367 j++; <<
368 } <<
369 } <<
370 <<
371 out.precision(10); <<
372 out.width(15); <<
373 out.setf(std::ofstream::left); <<
374 out << -1.f << ' '; <<
375 <<
376 out.precision(10); <<
377 out.width(15); <<
378 out.setf(std::ofstream::left); <<
379 out << -1.f << std::endl; <<
380 <<
381 k++; <<
382 } <<
383 <<
384 out.precision(10); <<
385 out.width(15); <<
386 out.setf(std::ofstream::left); <<
387 out << -2.f << ' '; <<
388 <<
389 out.precision(10); <<
390 out.width(15); <<
391 out.setf(std::ofstream::left); <<
392 out << -2.f << std::endl; <<
393 <<
394 return true; <<
395 } 201 }
396 202
397 //....oooOO0OOooo........oooOO0OOooo........oo << 203 const G4DataVector& G4ShellEMDataSet::GetData(G4int i) const
398 <<
399 void G4ShellEMDataSet::CleanUpComponents() <<
400 { 204 {
401 while (!components.empty()) << 205 const G4VEMDataSet* component = GetComponent(i);
402 { << 206 return (component->GetData(i));
403 if (components.back()) delete components <<
404 components.pop_back(); <<
405 } <<
406 } <<
407 <<
408 //....oooOO0OOooo........oooOO0OOooo........oo <<
409 <<
410 G4String G4ShellEMDataSet::FullFileName(const <<
411 { <<
412 const char* path = G4FindDataDir("G4LEDATA") <<
413 <<
414 if (!path) <<
415 { <<
416 G4Exception("G4ShellEMDataSet::FullFileNam <<
417 return ""; <<
418 } <<
419 <<
420 std::ostringstream fullFileName; <<
421 <<
422 fullFileName << path << '/' << fileName << z <<
423 <<
424 return G4String(fullFileName.str().c_str()); <<
425 } 207 }
426 208