Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4ENDFTapeRead.cc
Go to the documentation of this file.
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 * File: G4ENDFTapeRead.cc
28 * Author: B. Wendt (wendbryc@isu.edu)
29 *
30 * Created on September 6, 2011, 10:01 AM
31 */
32
33#include <fstream>
34#include <map>
35#include <vector>
36
37#include "globals.hh"
38#include "G4ParticleHPManager.hh"
39
40#include "G4ENDFTapeRead.hh"
41#include "G4ENDFYieldDataContainer.hh"
42#include "G4FFGDebuggingMacros.hh"
43#include "G4FFGDefaultValues.hh"
44#include "G4FFGEnumerations.hh"
45#include "G4TableTemplate.hh"
46
47G4ENDFTapeRead::
48G4ENDFTapeRead( G4String FileLocation,
49 G4String FileName,
50 G4FFGEnumerations::YieldType WhichYield,
51 G4FFGEnumerations::FissionCause /*WhichCause*/ )
52: /* Cause_(WhichCause),*/
53 Verbosity_(G4FFGDefaultValues::Verbosity),
54 YieldType_(WhichYield)
55{
56 // Initialize the class
57 Initialize(FileLocation + FileName);
58}
59
60G4ENDFTapeRead::
61G4ENDFTapeRead( G4String FileLocation,
62 G4String FileName,
63 G4FFGEnumerations::YieldType WhichYield,
64 G4FFGEnumerations::FissionCause /*WhichCause*/,
65 G4int Verbosity )
66: /*Cause_(WhichCause),*/
67 Verbosity_(Verbosity),
68 YieldType_(WhichYield)
69{
70 // Initialize the class
71 Initialize(FileLocation + FileName);
72}
73
74G4ENDFTapeRead::
75G4ENDFTapeRead( std::istringstream& dataStream,
76 G4FFGEnumerations::YieldType WhichYield,
77 G4FFGEnumerations::FissionCause /*WhichCause*/,
78 G4int Verbosity )
79: /*Cause_(WhichCause),*/
80 Verbosity_(Verbosity),
81 YieldType_(WhichYield)
82{
83 // Initialize the class
84 Initialize(dataStream);
85}
86
87void G4ENDFTapeRead::
88Initialize( G4String dataFile )
89{
90 std::istringstream dataStream(std::ios::in);
91 G4ParticleHPManager::GetInstance()->GetDataStream(dataFile, dataStream);
92
93 Initialize(dataStream);
94}
95
96void G4ENDFTapeRead::
97Initialize( std::istringstream& dataStream )
98{
99G4FFG_FUNCTIONENTER__
100
101 EnergyGroups_ = 0;
102 EnergyGroupValues_ = NULL;
103
104 YieldContainerTable_ = new G4TableTemplate< G4ENDFYieldDataContainer >;
105
106 try
107 {
108 ReadInData(dataStream);
109 } catch (std::exception & e)
110 {
111 delete YieldContainerTable_;
112
113 G4FFG_FUNCTIONLEAVE__
114 throw e;
115 }
116
117G4FFG_FUNCTIONLEAVE__
118}
119
120G4double* G4ENDFTapeRead::
121G4GetEnergyGroupValues( void )
122{
123G4FFG_FUNCTIONENTER__
124
125G4FFG_FUNCTIONLEAVE__
126 return EnergyGroupValues_;
127}
128
129G4int G4ENDFTapeRead::
130G4GetNumberOfEnergyGroups( void )
131{
132G4FFG_FUNCTIONENTER__
133
134G4FFG_FUNCTIONLEAVE__
135 return EnergyGroups_;
136}
137
138G4int G4ENDFTapeRead::
139G4GetNumberOfFissionProducts( void )
140{
141G4FFG_FUNCTIONENTER__
142
143 G4int NumberOfElements = (G4int)YieldContainerTable_->G4GetNumberOfElements();
144
145G4FFG_FUNCTIONLEAVE__
146 return NumberOfElements;
147}
148
149G4ENDFYieldDataContainer* G4ENDFTapeRead::
150G4GetYield( G4int WhichYield )
151{
152G4FFG_DATA_FUNCTIONENTER__
153
154 G4ENDFYieldDataContainer* Container = NULL;
155 if(WhichYield >= 0 && WhichYield < YieldContainerTable_->G4GetNumberOfElements())
156 {
157 Container = YieldContainerTable_->G4GetContainer(WhichYield);
158 }
159
160G4FFG_DATA_FUNCTIONLEAVE__
161 return Container;
162}
163
164void G4ENDFTapeRead::
165G4SetVerbosity(G4int WhatVerbosity)
166{
167G4FFG_FUNCTIONENTER__
168
169 this->Verbosity_ = WhatVerbosity;
170
171G4FFG_FUNCTIONLEAVE__
172}
173
174void G4ENDFTapeRead::
175ReadInData( std::istringstream& dataStream )
176{
177G4FFG_FUNCTIONENTER__
178
179 // Check if the file exists
180 if(!dataStream.good())
181 {
182 G4Exception("G4ENDFTapeRead::ReadInData()",
183 "Illegal file name",
184 JustWarning,
185 "Fission product data not available");
186
187 // TODO create/use a specialized exception
188 G4FFG_FUNCTIONLEAVE__
189 throw std::exception();
190 }
191
192 // Code to read in from a pure ENDF data tape.
193 // Commented out while pre-formatted Geant4 ENDF data is being used
194// G4int CurrentEnergyGroup = -1;
195// std::vector< G4double > NewDoubleVector;
196// std::vector< G4double > EnergyPoints;
197// std::vector< G4int > Product;
198// std::vector< G4FFGEnumerations::MetaState > MetaState;
199// std::vector< std::vector< G4double > > Yield;
200// std::vector< std::vector< G4double > > Error;
201// G4String DataBlock;
202// std::size_t InsertExponent;
203// G4double Parts[6];
204// G4double dummy;
205// G4int MAT;
206// G4int MF;
207// G4int MT;
208// G4int LN;
209// G4int Block;
210// G4int EmptyProduct;
211// G4int Location;
212// G4int ItemCounter = 0;
213// G4int FirstLineInEnergyGroup = 0;
214// G4int LastLineInEnergyGroup = 0;
215// G4bool FoundEnergyGroup = false;
216// G4bool FoundPID = false;
217//
218// while(getline(DataFile, Temp))
219// {
220// // Format the string so that it can be interpreted correctly
221// DataBlock = Temp.substr(0, 66);
222// Temp = Temp.substr(66);
223// InsertExponent = 0;
224// while((InsertExponent = DataBlock.find_first_of("-+", InsertExponent)) != G4String::npos)
225// {
226// DataBlock.insert(InsertExponent, 1, 'e');
227// InsertExponent += 2;
228// }
229// sscanf(DataBlock.c_str(), "%11le %11le %11le %11le %11le %11le",
230// &Parts[0], &Parts[1], &Parts[2], &Parts[3], &Parts[4], &Parts[5]);
231// sscanf(Temp.substr(0, 4).c_str(), "%i", &MAT);
232// sscanf(Temp.substr(4, 2).c_str(), "%i", &MF);
233// sscanf(Temp.substr(6, 3).c_str(), "%i", &MT);
234// sscanf(Temp.substr(9).c_str(), "%i", &LN);
235//
236// if(MT == YieldType_)
237// {
238// if(LN == 1)
239// {
240// if(FoundPID != true)
241// {
242// // The first line of an ENDF section for MT = 454 or 459
243// // always contains the parent PID
244// // This section can potentially be expanded to check and
245// // verify that it is the correct nucleus
246// FoundPID = true;
247//
248// continue;
249// }
250// } else if(FoundPID == true && FoundEnergyGroup == false)
251// {
252// // Skip this line if it is not the energy definition line
253// if(Parts[1] != 0 || Parts[3] != 0)
254// {
255// continue;
256// }
257//
258// // The first block is the incident neutron energy
259// // information.
260// // Check to make sure that it is spontaneous or neutron
261// // induced.
262// if(Cause_ == G4FFGEnumerations::NEUTRON_INDUCED)
263// {
264// if(Parts[0] != 0)
265// {
266// FoundEnergyGroup = true;
267// }
268// } else if(Cause_ == G4FFGEnumerations::SPONTANEOUS)
269// {
270// if(Parts[0] == 0)
271// {
272// FoundEnergyGroup = true;
273// }
274// } else
275// { // Maybe more fission causes here if added later
276// FoundEnergyGroup = false;
277// }
278//
279// if(FoundEnergyGroup == true)
280// {
281// // Convert to eV
282// Parts[0] *= eV;
283//
284// // Calculate the parameters
285// FirstLineInEnergyGroup = LN;
286// LastLineInEnergyGroup = FirstLineInEnergyGroup +
287// ceil(Parts[4] / 6.0);
288// ItemCounter = 0;
289// EmptyProduct = 0;
290//
291// // Initialize the data storage
292// CurrentEnergyGroup++;
293// EnergyPoints.push_back(Parts[0]);
294// Yield.push_back(NewDoubleVector);
295// Yield.back().resize(Product.size(), 0);
296// Error.push_back(NewDoubleVector);
297// Error.back().resize(Product.size(), 0);
298//
299// continue;
300// }
301// }
302//
303// if(LN > FirstLineInEnergyGroup && LN <= LastLineInEnergyGroup)
304// {
305// for(Block = 0; Block < 6; Block++)
306// {
307// if(EmptyProduct > 0)
308// {
309// EmptyProduct--;
310//
311// continue;
312// }
313// switch(ItemCounter % 4)
314// {
315// case 0:
316// // Determine if the block is empty
317// if(Parts[Block] == 0)
318// {
319// EmptyProduct = 3;
320//
321// continue;
322// }
323//
324// // Determine if this product is already loaded
325// for(Location = 0; Location < (signed)Product.size(); Location++)
326// {
327// if(Parts[Block] == Product.at(Location) &&
328// Parts[Block + 1] == MetaState.at(Location))
329// {
330// break;
331// }
332// }
333//
334// // The product hasn't been created yet
335// // Add it and initialize the other vectors
336// if(Location == (signed)Product.size())
337// {
338// Product.push_back(Parts[Block]);
339// MetaState.push_back((G4FFGEnumerations::MetaState)Parts[Block + 1]);
340// Yield.at(CurrentEnergyGroup).push_back(0.0);
341// Error.at(CurrentEnergyGroup).push_back(0.0);
342// }
343// break;
344//
345// case 2:
346// Yield.at(CurrentEnergyGroup).at(Location) = Parts[Block];
347// break;
348//
349// case 3:
350// Error.at(CurrentEnergyGroup).at(Location) = Parts[Block];
351// break;
352// }
353//
354// ItemCounter++;
355// }
356// }
357//
358// if (LN == LastLineInEnergyGroup)
359// {
360// FoundEnergyGroup = false;
361// }
362// }
363// }
364//
365// G4ENDFYieldDataContainer* NewDataContainer;
366// EnergyGroups_ = EnergyPoints.size();
367// EnergyGroupValues_ = new G4double[EnergyGroups_];
368// G4int NewProduct;
369// G4FFGEnumerations::MetaState NewMetaState;
370// G4double* NewYield = new G4double[EnergyGroups_];
371// G4double* NewError = new G4double[EnergyGroups_];
372//
373// for(G4int i = 0; i < EnergyGroups_; i++)
374// {
375// // Load the energy values
376// EnergyGroupValues_[i] = EnergyPoints.at(i);
377//
378// // Make all the vectors the same size
379// Yield[i].resize(maxSize, 0.0);
380// Error[i].resize(maxSize, 0.0);
381// }
382//
383// // Load the data into the yield table
384// for(ItemCounter = 0; ItemCounter < (signed)Product.size(); ItemCounter++)
385// {
386// NewProduct = Product.at(ItemCounter);
387// NewMetaState = MetaState.at(ItemCounter);
388//
389// for(CurrentEnergyGroup = 0; CurrentEnergyGroup < EnergyGroups_; CurrentEnergyGroup++)
390// {
391// NewYield[CurrentEnergyGroup] = Yield.at(CurrentEnergyGroup).at(ItemCounter);
392// NewYield[CurrentEnergyGroup] = Error.at(CurrentEnergyGroup).at(ItemCounter);
393// }
394//
395// NewDataContainer = YieldContainerTable_->G4GetNewContainer(EnergyGroups_ + 1);
396// NewDataContainer->SetProduct(NewProduct);
397// NewDataContainer->SetMetaState(NewMetaState);
398// NewDataContainer->SetYieldProbability(NewYield);
399// NewDataContainer->SetYieldError(NewError);
400// }
401//
402// delete[] NewYield;
403// delete[] NewError;
404
405 G4int MT;
406 G4bool correctMT;
407 G4int MF;
408 G4double dummy;
409 G4int blockCount;
410 G4int currentEnergy = 0;
411 G4double incidentEnergy;
412 G4int itemCount;
413 // TODO correctly implement the interpolation in the fission product yield
414 G4int interpolation;
415 G4int isotope;
416 G4int metastate;
417 G4int identifier;
418 G4double yield;
419 // "error" is included here in the event that errors are included in the future
420 G4double error = 0.0;
421 G4int maxSize = 0;
422
423 std::vector< G4double > projectileEnergies;
424 std::map< const G4int, std::pair< std::vector< G4double >, std::vector< G4double > > > intermediateData;
425 std::map< const G4int, std::pair< std::vector< G4double >, std::vector< G4double > > >::iterator dataIterator;
426
427 while(dataStream.good()) // Loop checking, 11.05.2015, T. Koi
428 {
429 dataStream >> MT >> MF >> dummy >> blockCount;
430
431 correctMT = MT == YieldType_;
432
433 for(G4int b = 0; b < blockCount; ++b)
434 {
435 dataStream >> incidentEnergy >> itemCount >> interpolation;
436 maxSize = maxSize >= itemCount ? maxSize : itemCount;
437
438 if(correctMT)
439 {
440 // Load in the energy of the projectile
441 projectileEnergies.push_back(incidentEnergy);
442 currentEnergy = G4int(projectileEnergies.size() - 1);
443 } else
444 {
445 // !!! Do not break since we need to parse through the !!!
446 // !!! entire data file for all possible energies !!!
447 }
448
449 for(G4int i = 0; i < itemCount; ++i)
450 {
451 dataStream >> isotope >> metastate >> yield;
452
453 if(correctMT)
454 {
455 identifier = isotope * 10 + metastate;
456
457 dataIterator = intermediateData.insert(std::make_pair(
458 identifier,
459 std::make_pair(
460 std::vector< G4double >(projectileEnergies.size(), 0.0),
461 std::vector< G4double >(projectileEnergies.size(), 0.0)))).first;
462
463 if(dataIterator->second.first.size() < projectileEnergies.size())
464 {
465 dataIterator->second.first.resize(projectileEnergies.size());
466 dataIterator->second.second.resize(projectileEnergies.size());
467 }
468
469 dataIterator->second.first[currentEnergy] = yield;
470 dataIterator->second.second[currentEnergy] = error;
471 } else
472 {
473 // !!! Do not break since we need to parse through the !!!
474 // !!! entire data file for all possible energies !!!
475 }
476 }
477 }
478 }
479
480 G4ENDFYieldDataContainer* NewDataContainer;
481 EnergyGroups_ = (G4int)projectileEnergies.size();
482 EnergyGroupValues_ = new G4double[EnergyGroups_];
483 G4int NewProduct;
484 G4FFGEnumerations::MetaState NewMetaState;
485 G4double* NewYield = new G4double[EnergyGroups_];
486 G4double* NewError = new G4double[EnergyGroups_];
487
488 for(G4int energyGroup = 0; energyGroup < EnergyGroups_; energyGroup++)
489 {
490 // Load the energy values
491 EnergyGroupValues_[energyGroup] = projectileEnergies[energyGroup];
492 }
493
494 // Load the data into the yield table
495 for(dataIterator = intermediateData.begin(); dataIterator != intermediateData.end(); ++dataIterator)
496 {
497 identifier = dataIterator->first;
498 metastate = identifier % 10;
499 switch(metastate)
500 {
501 case 1:
502 NewMetaState = G4FFGEnumerations::META_1;
503 break;
504
505 case 2:
506 NewMetaState = G4FFGEnumerations::META_2;
507 break;
508
509 default:
510 G4Exception("G4ENDFTapeRead::ReadInData()",
511 "Unsupported state",
512 JustWarning,
513 "Unsupported metastable state supplied in fission yield data. Defaulting to the ground state");
514 // Fall through
515 case 0:
516 NewMetaState = G4FFGEnumerations::GROUND_STATE;
517 break;
518 }
519 NewProduct = (identifier - metastate) / 10;
520
521 for(G4int energyGroup = 0; energyGroup < EnergyGroups_; energyGroup++)
522 {
523 if(energyGroup < (signed)dataIterator->second.first.size())
524 {
525 yield = dataIterator->second.first[energyGroup];
526 error = dataIterator->second.second[energyGroup];
527 } else
528 {
529 yield = 0.0;
530 error = 0.0;
531 }
532
533 NewYield[energyGroup] = yield;
534 NewError[energyGroup] = error;
535 }
536
537 NewDataContainer = YieldContainerTable_->G4GetNewContainer(EnergyGroups_);
538 NewDataContainer->SetProduct(NewProduct);
539 NewDataContainer->SetMetaState(NewMetaState);
540 NewDataContainer->SetYieldProbability(NewYield);
541 NewDataContainer->SetYieldError(NewError);
542 }
543
544 delete[] NewYield;
545 delete[] NewError;
546
547G4FFG_FUNCTIONLEAVE__
548}
549
550G4ENDFTapeRead::
551~G4ENDFTapeRead( void )
552{
553G4FFG_FUNCTIONENTER__
554
555 delete[] EnergyGroupValues_;
556 delete YieldContainerTable_;
557
558G4FFG_FUNCTIONLEAVE__
559}
560
JustWarning
@ JustWarning
Definition: G4ExceptionSeverity.hh:65
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:59
G4FFG_DATA_FUNCTIONENTER__
#define G4FFG_DATA_FUNCTIONENTER__
Definition: G4FFGDebuggingMacros.hh:218
G4FFG_DATA_FUNCTIONLEAVE__
#define G4FFG_DATA_FUNCTIONLEAVE__
Definition: G4FFGDebuggingMacros.hh:222
G4FFG_FUNCTIONLEAVE__
#define G4FFG_FUNCTIONLEAVE__
Definition: G4FFGDebuggingMacros.hh:219
G4FFG_FUNCTIONENTER__
#define G4FFG_FUNCTIONENTER__
Definition: G4FFGDebuggingMacros.hh:215
G4double
double G4double
Definition: G4Types.hh:83
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
G4ENDFTapeRead::G4GetNumberOfFissionProducts
G4int G4GetNumberOfFissionProducts(void)
Definition: G4ENDFTapeRead.cc:139
G4ENDFTapeRead::Initialize
void Initialize(G4String dataFile)
Definition: G4ENDFTapeRead.cc:88
G4ENDFTapeRead::G4SetVerbosity
void G4SetVerbosity(G4int WhatVerbosity)
Definition: G4ENDFTapeRead.cc:165
G4ENDFTapeRead::G4GetEnergyGroupValues
G4double * G4GetEnergyGroupValues(void)
Definition: G4ENDFTapeRead.cc:121
G4ENDFTapeRead::G4GetYield
G4ENDFYieldDataContainer * G4GetYield(G4int WhichYield)
Definition: G4ENDFTapeRead.cc:150
G4ENDFTapeRead::G4ENDFTapeRead
G4ENDFTapeRead(G4String FileLocation, G4String FileName, G4FFGEnumerations::YieldType WhichYield, G4FFGEnumerations::FissionCause WhichCause)
Definition: G4ENDFTapeRead.cc:48
G4ENDFTapeRead::G4GetNumberOfEnergyGroups
G4int G4GetNumberOfEnergyGroups(void)
Definition: G4ENDFTapeRead.cc:130
G4ENDFYieldDataContainer::SetMetaState
void SetMetaState(G4FFGEnumerations::MetaState MetaState)
Definition: G4ENDFYieldDataContainer.cc:81
G4ENDFYieldDataContainer::SetYieldError
void SetYieldError(G4double *YieldError)
Definition: G4ENDFYieldDataContainer.cc:93
G4ENDFYieldDataContainer::SetYieldProbability
void SetYieldProbability(G4double *YieldProbability)
Definition: G4ENDFYieldDataContainer.cc:99
G4ParticleHPManager::GetInstance
static G4ParticleHPManager * GetInstance()
Definition: G4ParticleHPManager.cc:73
G4ParticleHPManager::GetDataStream
void GetDataStream(G4String, std::istringstream &iss)
Definition: G4ParticleHPManager.cc:98
G4TableTemplate::G4GetNewContainer
T * G4GetNewContainer(void)
Definition: G4TableTemplate.hh:106
G4TableTemplate::G4GetNumberOfElements
G4long G4GetNumberOfElements(void)
Definition: G4TableTemplate.hh:124
G4TableTemplate::G4GetContainer
T * G4GetContainer(unsigned int WhichContainer)
Definition: G4TableTemplate.hh:94