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