| 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 * 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&
48 G4FFGEnumeratio
49 G4FFGEnumeratio
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&
59 G4FFGEnumeratio
60 G4FFGEnumeratio
61 : /*Cause_(WhichCause),*/
62 Verbosity_(Verbosity),
63 YieldType_(WhichYield)
64 {
65 // Initialize the class
66 Initialize(FileLocation + FileName);
67 }
68
69 G4ENDFTapeRead::G4ENDFTapeRead(std::istringstr
70 G4FFGEnumeratio
71 G4FFGEnumeratio
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
81 {
82 std::istringstream dataStream(std::ios::in);
83 G4ParticleHPManager::GetInstance()->GetDataS
84
85 Initialize(dataStream);
86 }
87
88 void G4ENDFTapeRead::Initialize(std::istringst
89 {
90 G4FFG_FUNCTIONENTER__
91
92 EnergyGroups_ = 0;
93 EnergyGroupValues_ = nullptr;
94
95 YieldContainerTable_ = new G4TableTemplate<G
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::G4GetEnergyGroupValu
111 {
112 G4FFG_FUNCTIONENTER__
113
114 G4FFG_FUNCTIONLEAVE__
115 return EnergyGroupValues_;
116 }
117
118 G4int G4ENDFTapeRead::G4GetNumberOfEnergyGroup
119 {
120 G4FFG_FUNCTIONENTER__
121
122 G4FFG_FUNCTIONLEAVE__
123 return EnergyGroups_;
124 }
125
126 G4int G4ENDFTapeRead::G4GetNumberOfFissionProd
127 {
128 G4FFG_FUNCTIONENTER__
129
130 auto NumberOfElements = (G4int)YieldContaine
131
132 G4FFG_FUNCTIONLEAVE__
133 return NumberOfElements;
134 }
135
136 G4ENDFYieldDataContainer* G4ENDFTapeRead::G4Ge
137 {
138 G4FFG_DATA_FUNCTIONENTER__
139
140 G4ENDFYieldDataContainer* Container = nullpt
141 if (WhichYield >= 0 && WhichYield < YieldCon
142 Container = YieldContainerTable_->G4GetCon
143 }
144
145 G4FFG_DATA_FUNCTIONLEAVE__
146 return Container;
147 }
148
149 void G4ENDFTapeRead::G4SetVerbosity(G4int What
150 {
151 G4FFG_FUNCTIONENTER__
152
153 this->Verbosity_ = WhatVerbosity;
154
155 G4FFG_FUNCTIONLEAVE__
156 }
157
158 void G4ENDFTapeRead::ReadInData(std::istringst
159 {
160 G4FFG_FUNCTIONENTER__
161
162 // Check if the file exists
163 if (!dataStream.good()) {
164 G4Exception("G4ENDFTapeRead::ReadInData()"
165 "Fission product data not avai
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 tap
173 // Commented out while pre-formatted Geant4
174 // G4int CurrentEnergyGroup = -1;
175 // std::vector< G4double > NewDoubleVecto
176 // std::vector< G4double > EnergyPoints;
177 // std::vector< G4int > Product;
178 // std::vector< G4FFGEnumerations::MetaSt
179 // std::vector< std::vector< G4double > >
180 // std::vector< std::vector< G4double > >
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 ca
201 // DataBlock = Temp.substr(0, 66);
202 // Temp = Temp.substr(66);
203 // InsertExponent = 0;
204 // while((InsertExponent = DataBlock.
205 // G4String::npos)
206 // {
207 // DataBlock.insert(InsertExponen
208 // InsertExponent += 2;
209 // }
210 // sscanf(DataBlock.c_str(), "%11le %
211 // &Parts[0], &Parts[1], &Parts[2
212 // sscanf(Temp.substr(0, 4).c_str(),
213 // sscanf(Temp.substr(4, 2).c_str(),
214 // sscanf(Temp.substr(6, 3).c_str(),
215 // sscanf(Temp.substr(9).c_str(), "%i
216 //
217 // if(MT == YieldType_)
218 // {
219 // if(LN == 1)
220 // {
221 // if(FoundPID != true)
222 // {
223 // // The first line of a
224 // // always contains the
225 // // This section can po
226 // // verify that it is t
227 // FoundPID = true;
228 //
229 // continue;
230 // }
231 // } else if(FoundPID == true &&
232 // {
233 // // Skip this line if it is
234 // if(Parts[1] != 0 || Parts[
235 // {
236 // continue;
237 // }
238 //
239 // // The first block is the
240 // // information.
241 // // Check to make sure that
242 // // induced.
243 // if(Cause_ == G4FFGEnumerat
244 // {
245 // if(Parts[0] != 0)
246 // {
247 // FoundEnergyGroup =
248 // }
249 // } else if(Cause_ == G4FFGE
250 // {
251 // if(Parts[0] == 0)
252 // {
253 // FoundEnergyGroup =
254 // }
255 // } else
256 // { // Maybe more fission ca
257 // FoundEnergyGroup = fal
258 // }
259 //
260 // if(FoundEnergyGroup == tru
261 // {
262 // // Convert to eV
263 // Parts[0] *= eV;
264 //
265 // // Calculate the param
266 // FirstLineInEnergyGroup
267 // LastLineInEnergyGroup
268 // ceil(Parts[4] / 6.
269 // ItemCounter = 0;
270 // EmptyProduct = 0;
271 //
272 // // Initialize the data
273 // CurrentEnergyGroup++;
274 // EnergyPoints.push_back
275 // Yield.push_back(NewDou
276 // Yield.back().resize(Pr
277 // Error.push_back(NewDou
278 // Error.back().resize(Pr
279 //
280 // continue;
281 // }
282 // }
283 //
284 // if(LN > FirstLineInEnergyGroup
285 // {
286 // for(Block = 0; Block < 6;
287 // {
288 // if(EmptyProduct > 0)
289 // {
290 // EmptyProduct--;
291 //
292 // continue;
293 // }
294 // switch(ItemCounter % 4
295 // {
296 // case 0:
297 // // Determine i
298 // if(Parts[Block
299 // {
300 // EmptyProdu
301 //
302 // continue;
303 // }
304 //
305 // // Determine i
306 // for(Location =
307 // {
308 // if(Parts[B
309 // Parts[B
310 // {
311 // break;
312 // }
313 // }
314 //
315 // // The product
316 // // Add it and
317 // if(Location ==
318 // {
319 // Product.pu
320 // MetaState.
321 // 1]); Yield
322 // Error.at(C
323 // }
324 // break;
325 //
326 // case 2:
327 // Yield.at(Curre
328 // break;
329 //
330 // case 3:
331 // Error.at(Curre
332 // break;
333 // }
334 //
335 // ItemCounter++;
336 // }
337 // }
338 //
339 // if (LN == LastLineInEnergyGrou
340 // {
341 // FoundEnergyGroup = false;
342 // }
343 // }
344 // }
345 //
346 // G4ENDFYieldDataContainer* NewDataConta
347 // EnergyGroups_ = EnergyPoints.size();
348 // EnergyGroupValues_ = new G4double[Ener
349 // G4int NewProduct;
350 // G4FFGEnumerations::MetaState NewMetaSt
351 // G4double* NewYield = new G4double[Ener
352 // G4double* NewError = new G4double[Ener
353 //
354 // for(G4int i = 0; i < EnergyGroups_; i+
355 // {
356 // // Load the energy values
357 // EnergyGroupValues_[i] = EnergyPoin
358 //
359 // // Make all the vectors the same s
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 < (si
366 // {
367 // NewProduct = Product.at(ItemCounte
368 // NewMetaState = MetaState.at(ItemCo
369 //
370 // for(CurrentEnergyGroup = 0; Curren
371 // {
372 // NewYield[CurrentEnergyGroup] =
373 // NewYield[CurrentEnergyGroup] =
374 // }
375 //
376 // NewDataContainer = YieldContainerT
377 // NewDataContainer->SetProduct(NewPr
378 // NewDataContainer->SetMetaState(New
379 // NewDataContainer->SetYieldProbabil
380 // NewDataContainer->SetYieldError(Ne
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 interpolatio
395 G4int interpolation;
396 G4int isotope;
397 G4int metastate;
398 G4int identifier;
399 G4double yield;
400 // "error" is included here in the event tha
401 G4double error = 0.0;
402 G4int maxSize = 0;
403
404 std::vector<G4double> projectileEnergies;
405 std::map<const G4int, std::pair<std::vector<
406 std::map<const G4int, std::pair<std::vector<
407 dataIterator;
408
409 while (dataStream.good()) // Loop checking,
410 {
411 dataStream >> MT >> MF >> dummy >> blockCo
412
413 correctMT = MT == YieldType_;
414
415 for (G4int b = 0; b < blockCount; ++b) {
416 dataStream >> incidentEnergy >> itemCoun
417 maxSize = maxSize >= itemCount ? maxSize
418
419 if (correctMT) {
420 // Load in the energy of the projectil
421 projectileEnergies.push_back(incidentE
422 currentEnergy = G4int(projectileEnergi
423 }
424 else {
425 // !!! Do not break since we need to p
426 // !!! entire data file for all possib
427 }
428
429 for (G4int i = 0; i < itemCount; ++i) {
430 dataStream >> isotope >> metastate >>
431
432 if (correctMT) {
433 identifier = isotope * 10 + metastat
434
435 dataIterator =
436 intermediateData
437 .insert(std::make_pair(
438 identifier, std::make_pair(std
439 std
440 .first;
441
442 if (dataIterator->second.first.size(
443 dataIterator->second.first.resize(
444 dataIterator->second.second.resize
445 }
446
447 dataIterator->second.first[currentEn
448 dataIterator->second.second[currentE
449 }
450 else {
451 // !!! Do not break since we need to
452 // !!! entire data file for all poss
453 }
454 }
455 }
456 }
457
458 G4ENDFYieldDataContainer* NewDataContainer;
459 EnergyGroups_ = (G4int)projectileEnergies.si
460 EnergyGroupValues_ = new G4double[EnergyGrou
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 < En
467 // Load the energy values
468 EnergyGroupValues_[energyGroup] = projecti
469 }
470
471 // Load the data into the yield table
472 for (dataIterator = intermediateData.begin()
473 ++dataIterator)
474 {
475 identifier = dataIterator->first;
476 metastate = identifier % 10;
477 switch (metastate) {
478 case 1:
479 NewMetaState = G4FFGEnumerations::META
480 break;
481
482 case 2:
483 NewMetaState = G4FFGEnumerations::META
484 break;
485
486 default:
487 G4Exception("G4ENDFTapeRead::ReadInDat
488 "Unsupported metastable st
489 "the ground state");
490 // Fall through
491 case 0:
492 NewMetaState = G4FFGEnumerations::GROU
493 break;
494 }
495 NewProduct = (identifier - metastate) / 10
496
497 for (G4int energyGroup = 0; energyGroup <
498 if (energyGroup < (signed)dataIterator->
499 yield = dataIterator->second.first[ene
500 error = dataIterator->second.second[en
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_->G
512 NewDataContainer->SetProduct(NewProduct);
513 NewDataContainer->SetMetaState(NewMetaStat
514 NewDataContainer->SetYieldProbability(NewY
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