source: trunk/source/processes/electromagnetic/lowenergy/src/G4DNACrossSectionDataSet.cc

Last change on this file was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 14.5 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
28// $Id: G4DNACrossSectionDataSet.cc,v 1.11 2009/11/12 10:05:30 sincerti Exp $
29// GEANT4 tag $Name: geant4-09-04-ref-00 $
30//
31// Author: Riccardo Capra <capra@ge.infn.it>
32// Code review by MGP October 2007: removed inheritance from concrete class
33// more improvements needed
34//
35// History:
36// -----------
37// 30 Jun 2005 RC Created
38// 14 Oct 2007 MGP Removed inheritance from concrete class G4ShellEMDataSet
39//
40// 15 Jul 2009 N.A.Karakatsanis
41//
42// - LoadNonLogData method was created to load only the non-logarithmic data from G4EMLOW
43// dataset. It is essentially performing the data loading operations as in the past.
44//
45// - LoadData method was revised in order to calculate the logarithmic values of the data
46// It retrieves the data values from the G4EMLOW data files but, then, calculates the
47// respective log values and loads them to seperate data structures.
48//
49// - SetLogEnergiesData method was cretaed to set logarithmic values to G4 data vectors.
50// The EM data sets, initialized this way, contain both non-log and log values.
51// These initialized data sets can enhance the computing performance of data interpolation
52// operations
53//
54//
55// -------------------------------------------------------------------
56
57
58#include "G4DNACrossSectionDataSet.hh"
59#include "G4VDataSetAlgorithm.hh"
60#include "G4EMDataSet.hh"
61#include <vector>
62#include <fstream>
63#include <sstream>
64
65
66G4DNACrossSectionDataSet::G4DNACrossSectionDataSet(G4VDataSetAlgorithm* argAlgorithm,
67 G4double argUnitEnergies,
68 G4double argUnitData)
69 :
70 algorithm(argAlgorithm), unitEnergies(argUnitEnergies), unitData(argUnitData)
71{
72 z = 0;
73
74}
75
76
77G4DNACrossSectionDataSet::~G4DNACrossSectionDataSet()
78{
79 CleanUpComponents();
80
81 if (algorithm)
82 delete algorithm;
83}
84
85G4bool G4DNACrossSectionDataSet::LoadData(const G4String & argFileName)
86{
87 CleanUpComponents();
88
89 G4String fullFileName(FullFileName(argFileName));
90 std::ifstream in(fullFileName, std::ifstream::binary|std::ifstream::in);
91
92 if (!in.is_open())
93 {
94 G4String message("G4DNACrossSectionDataSet::LoadData - data file \"");
95 message+=fullFileName;
96 message+="\" not found";
97 G4Exception(message);
98 return false;
99 }
100
101 std::vector<G4DataVector *> columns;
102 std::vector<G4DataVector *> log_columns;
103
104 std::stringstream *stream(new std::stringstream);
105 char c;
106 G4bool comment(false);
107 G4bool space(true);
108 G4bool first(true);
109
110 try
111 {
112 while (!in.eof())
113 {
114 in.get(c);
115
116 switch (c)
117 {
118 case '\r':
119 case '\n':
120 if (!first)
121 {
122 unsigned long i(0);
123 G4double value;
124
125 while (!stream->eof())
126 {
127 (*stream) >> value;
128
129 while (i>=columns.size())
130 {
131 columns.push_back(new G4DataVector);
132 log_columns.push_back(new G4DataVector);
133 }
134
135 columns[i]->push_back(value);
136
137// N. A. Karakatsanis
138// A condition is applied to check if negative or zero values are present in the dataset.
139// If yes, then a near-zero value is applied to allow the computation of the logarithmic value
140// If a value is zero, this simplification is acceptable
141// If a value is negative, then it is not acceptable and the data of the particular column of
142// logarithmic values should not be used by interpolation methods.
143//
144// Therefore, G4LogLogInterpolation and G4LinLogLogInterpolation should not be used if negative values are present.
145// Instead, G4LinInterpolation is safe in every case
146// SemiLogInterpolation is safe only if the energy columns are non-negative
147// G4LinLogInterpolation is safe only if the cross section data columns are non-negative
148
149 if (value <=0.) value = 1e-300;
150 log_columns[i]->push_back(std::log10(value));
151
152 i++;
153 }
154
155 delete stream;
156 stream=new std::stringstream;
157 }
158
159 first=true;
160 comment=false;
161 space=true;
162 break;
163
164 case '#':
165 comment=true;
166 break;
167
168 case '\t':
169 c=' ';
170 case ' ':
171 if (space)
172 break;
173 default:
174 if (comment)
175 break;
176
177 if (c==' ')
178 space=true;
179 else
180 {
181 if (space && (!first))
182 (*stream) << ' ';
183
184 first=false;
185 (*stream) << c;
186 space=false;
187 }
188 }
189 }
190 }
191 catch(const std::ios::failure &e)
192 {
193 // some implementations of STL could throw a "failture" exception
194 // when read wants read characters after end of file
195 }
196
197 delete stream;
198
199 std::vector<G4DataVector *>::size_type maxI(columns.size());
200
201 if (maxI<2)
202 {
203 G4String message("G4DNACrossSectionDataSet::LoadData - data file \"");
204 message+=fullFileName;
205 message+="\" should have at least two columns";
206 G4Exception(message);
207 return false;
208 }
209
210 std::vector<G4DataVector*>::size_type i(1);
211 while (i<maxI)
212 {
213 G4DataVector::size_type maxJ(columns[i]->size());
214
215 if (maxJ!=columns[0]->size())
216 {
217 G4String message("G4DNACrossSectionDataSet::LoadData - data file \"");
218 message+=fullFileName;
219 message+="\" has lines with a different number of columns";
220 G4Exception(message);
221 return false;
222 }
223
224 G4DataVector::size_type j(0);
225
226 G4DataVector *argEnergies=new G4DataVector;
227 G4DataVector *argData=new G4DataVector;
228 G4DataVector *argLogEnergies=new G4DataVector;
229 G4DataVector *argLogData=new G4DataVector;
230
231 while(j<maxJ)
232 {
233 argEnergies->push_back(columns[0]->operator[] (j)*GetUnitEnergies());
234 argData->push_back(columns[i]->operator[] (j)*GetUnitData());
235 argLogEnergies->push_back(log_columns[0]->operator[] (j) + std::log10(GetUnitEnergies()));
236 argLogData->push_back(log_columns[i]->operator[] (j) + std::log10(GetUnitData()));
237 j++;
238 }
239
240 AddComponent(new G4EMDataSet(i-1, argEnergies, argData, argLogEnergies, argLogData, GetAlgorithm()->Clone(), GetUnitEnergies(), GetUnitData()));
241
242 i++;
243 }
244
245 i=maxI;
246 while (i>0)
247 {
248 i--;
249 delete columns[i];
250 delete log_columns[i];
251 }
252
253 return true;
254}
255
256
257G4bool G4DNACrossSectionDataSet::LoadNonLogData(const G4String & argFileName)
258{
259 CleanUpComponents();
260
261 G4String fullFileName(FullFileName(argFileName));
262 std::ifstream in(fullFileName, std::ifstream::binary|std::ifstream::in);
263
264 if (!in.is_open())
265 {
266 G4String message("G4DNACrossSectionDataSet::LoadData - data file \"");
267 message+=fullFileName;
268 message+="\" not found";
269 G4Exception(message);
270 return false;
271 }
272
273 std::vector<G4DataVector *> columns;
274
275 std::stringstream *stream(new std::stringstream);
276 char c;
277 G4bool comment(false);
278 G4bool space(true);
279 G4bool first(true);
280
281 try
282 {
283 while (!in.eof())
284 {
285 in.get(c);
286
287 switch (c)
288 {
289 case '\r':
290 case '\n':
291 if (!first)
292 {
293 unsigned long i(0);
294 G4double value;
295
296 while (!stream->eof())
297 {
298 (*stream) >> value;
299
300 while (i>=columns.size())
301 {
302 columns.push_back(new G4DataVector);
303 }
304
305 columns[i]->push_back(value);
306
307 i++;
308 }
309
310 delete stream;
311 stream=new std::stringstream;
312 }
313
314 first=true;
315 comment=false;
316 space=true;
317 break;
318
319 case '#':
320 comment=true;
321 break;
322
323 case '\t':
324 c=' ';
325 case ' ':
326 if (space)
327 break;
328 default:
329 if (comment)
330 break;
331
332 if (c==' ')
333 space=true;
334 else
335 {
336 if (space && (!first))
337 (*stream) << ' ';
338
339 first=false;
340 (*stream) << c;
341 space=false;
342 }
343 }
344 }
345 }
346 catch(const std::ios::failure &e)
347 {
348 // some implementations of STL could throw a "failture" exception
349 // when read wants read characters after end of file
350 }
351
352 delete stream;
353
354 std::vector<G4DataVector *>::size_type maxI(columns.size());
355
356 if (maxI<2)
357 {
358 G4String message("G4DNACrossSectionDataSet::LoadData - data file \"");
359 message+=fullFileName;
360 message+="\" should have at least two columns";
361 G4Exception(message);
362 return false;
363 }
364
365 std::vector<G4DataVector*>::size_type i(1);
366 while (i<maxI)
367 {
368 G4DataVector::size_type maxJ(columns[i]->size());
369
370 if (maxJ!=columns[0]->size())
371 {
372 G4String message("G4DNACrossSectionDataSet::LoadData - data file \"");
373 message+=fullFileName;
374 message+="\" has lines with a different number of columns";
375 G4Exception(message);
376 return false;
377 }
378
379 G4DataVector::size_type j(0);
380
381 G4DataVector *argEnergies=new G4DataVector;
382 G4DataVector *argData=new G4DataVector;
383
384 while(j<maxJ)
385 {
386 argEnergies->push_back(columns[0]->operator[] (j)*GetUnitEnergies());
387 argData->push_back(columns[i]->operator[] (j)*GetUnitData());
388 j++;
389 }
390
391 AddComponent(new G4EMDataSet(i-1, argEnergies, argData, GetAlgorithm()->Clone(), GetUnitEnergies(), GetUnitData()));
392
393 i++;
394 }
395
396 i=maxI;
397 while (i>0)
398 {
399 i--;
400 delete columns[i];
401 }
402
403 return true;
404}
405
406
407G4bool G4DNACrossSectionDataSet::SaveData(const G4String & argFileName) const
408{
409 const size_t n(NumberOfComponents());
410
411 if (n==0)
412 {
413 G4Exception("G4EMDataSet::SaveData - expected at least one component");
414 return false;
415 }
416
417 G4String fullFileName(FullFileName(argFileName));
418 std::ofstream out(fullFileName);
419
420 if (!out.is_open())
421 {
422 G4String message("G4EMDataSet::SaveData - cannot open \"");
423 message+=fullFileName;
424 message+="\"";
425 G4Exception(message);
426 return false;
427 }
428
429 G4DataVector::const_iterator iEnergies(GetComponent(0)->GetEnergies(0).begin());
430 G4DataVector::const_iterator iEnergiesEnd(GetComponent(0)->GetEnergies(0).end());
431 G4DataVector::const_iterator * iData(new G4DataVector::const_iterator[n]);
432
433 size_t k(n);
434
435 while (k>0)
436 {
437 k--;
438 iData[k]=GetComponent(k)->GetData(0).begin();
439 }
440
441 while (iEnergies!=iEnergiesEnd)
442 {
443 out.precision(10);
444 out.width(15);
445 out.setf(std::ofstream::left);
446 out << ((*iEnergies)/GetUnitEnergies());
447
448 k=0;
449
450 while (k<n)
451 {
452 out << ' ';
453 out.precision(10);
454 out.width(15);
455 out.setf(std::ofstream::left);
456 out << ((*(iData[k]))/GetUnitData());
457
458 iData[k]++;
459 k++;
460 }
461
462 out << std::endl;
463
464 iEnergies++;
465 }
466
467 delete[] iData;
468
469 return true;
470}
471
472
473G4String G4DNACrossSectionDataSet::FullFileName(const G4String& argFileName) const
474{
475 char* path = getenv("G4LEDATA");
476 if (!path)
477 G4Exception("G4DNACrossSectionDataSet::FullFileName - G4LEDATA environment variable not set");
478
479 std::ostringstream fullFileName;
480
481 fullFileName << path << "/" << argFileName << ".dat";
482
483 return G4String(fullFileName.str().c_str());
484}
485
486
487G4double G4DNACrossSectionDataSet::FindValue(G4double argEnergy, G4int /* argComponentId */) const
488{
489 // Returns the sum over the shells corresponding to e
490 G4double value = 0.;
491
492 std::vector<G4VEMDataSet *>::const_iterator i(components.begin());
493 std::vector<G4VEMDataSet *>::const_iterator end(components.end());
494
495 while (i!=end)
496 {
497 value+=(*i)->FindValue(argEnergy);
498 i++;
499 }
500
501 return value;
502}
503
504
505void G4DNACrossSectionDataSet::PrintData(void) const
506{
507 const size_t n(NumberOfComponents());
508
509 G4cout << "The data set has " << n << " components" << G4endl;
510 G4cout << G4endl;
511
512 size_t i(0);
513
514 while (i<n)
515 {
516 G4cout << "--- Component " << i << " ---" << G4endl;
517 GetComponent(i)->PrintData();
518 i++;
519 }
520}
521
522
523void G4DNACrossSectionDataSet::SetEnergiesData(G4DataVector* argEnergies,
524 G4DataVector* argData,
525 G4int argComponentId)
526{
527 G4VEMDataSet * component(components[argComponentId]);
528
529 if (component)
530 {
531 component->SetEnergiesData(argEnergies, argData, 0);
532 return;
533 }
534
535 std::ostringstream message;
536 message << "G4DNACrossSectionDataSet::SetEnergiesData - component " << argComponentId << " not found";
537
538 G4Exception(message.str().c_str());
539}
540
541
542void G4DNACrossSectionDataSet::SetLogEnergiesData(G4DataVector* argEnergies,
543 G4DataVector* argData,
544 G4DataVector* argLogEnergies,
545 G4DataVector* argLogData,
546 G4int argComponentId)
547{
548 G4VEMDataSet * component(components[argComponentId]);
549
550 if (component)
551 {
552 component->SetLogEnergiesData(argEnergies, argData, argLogEnergies, argLogData, 0);
553 return;
554 }
555
556 std::ostringstream message;
557 message << "G4DNACrossSectionDataSet::SetLogEnergiesData - component " << argComponentId << " not found";
558
559 G4Exception(message.str().c_str());
560}
561
562
563void G4DNACrossSectionDataSet::CleanUpComponents()
564{
565 while (!components.empty())
566 {
567 if (components.back()) delete components.back();
568 components.pop_back();
569 }
570}
571
572
Note: See TracBrowser for help on using the repository browser.