| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 /* 26 /*
27 * File: G4ENDFTapeRead.cc 27 * File: G4ENDFTapeRead.cc
28 * Author: B. Wendt (wendbryc@isu.edu) 28 * Author: B. Wendt (wendbryc@isu.edu)
29 * 29 *
30 * Created on September 6, 2011, 10:01 AM 30 * Created on September 6, 2011, 10:01 AM
31 */ 31 */
32 32
33 #include "G4ENDFTapeRead.hh" << 33 #include <fstream>
>> 34 #include <map>
>> 35 #include <vector>
>> 36
>> 37 #include "globals.hh"
>> 38 #include "G4ParticleHPManager.hh"
34 39
>> 40 #include "G4ENDFTapeRead.hh"
35 #include "G4ENDFYieldDataContainer.hh" 41 #include "G4ENDFYieldDataContainer.hh"
36 #include "G4FFGDebuggingMacros.hh" 42 #include "G4FFGDebuggingMacros.hh"
37 #include "G4FFGDefaultValues.hh" 43 #include "G4FFGDefaultValues.hh"
38 #include "G4FFGEnumerations.hh" 44 #include "G4FFGEnumerations.hh"
39 #include "G4ParticleHPManager.hh" <<
40 #include "G4TableTemplate.hh" 45 #include "G4TableTemplate.hh"
41 #include "globals.hh" <<
42 <<
43 #include <fstream> <<
44 #include <map> <<
45 #include <vector> <<
46 46
47 G4ENDFTapeRead::G4ENDFTapeRead(const G4String& << 47 G4ENDFTapeRead::
48 G4FFGEnumeratio << 48 G4ENDFTapeRead( G4String FileLocation,
49 G4FFGEnumeratio << 49 G4String FileName,
50 : /* Cause_(WhichCause),*/ << 50 G4FFGEnumerations::YieldType WhichYield,
>> 51 G4FFGEnumerations::FissionCause /*WhichCause*/ )
>> 52 : /* Cause_(WhichCause),*/
51 Verbosity_(G4FFGDefaultValues::Verbosity), 53 Verbosity_(G4FFGDefaultValues::Verbosity),
52 YieldType_(WhichYield) 54 YieldType_(WhichYield)
53 { 55 {
54 // Initialize the class << 56 // Initialize the class
55 Initialize(FileLocation + FileName); << 57 Initialize(FileLocation + FileName);
56 } 58 }
57 59
58 G4ENDFTapeRead::G4ENDFTapeRead(const G4String& << 60 G4ENDFTapeRead::
59 G4FFGEnumeratio << 61 G4ENDFTapeRead( G4String FileLocation,
60 G4FFGEnumeratio << 62 G4String FileName,
61 : /*Cause_(WhichCause),*/ << 63 G4FFGEnumerations::YieldType WhichYield,
>> 64 G4FFGEnumerations::FissionCause /*WhichCause*/,
>> 65 G4int Verbosity )
>> 66 : /*Cause_(WhichCause),*/
62 Verbosity_(Verbosity), 67 Verbosity_(Verbosity),
63 YieldType_(WhichYield) 68 YieldType_(WhichYield)
64 { 69 {
65 // Initialize the class << 70 // Initialize the class
66 Initialize(FileLocation + FileName); << 71 Initialize(FileLocation + FileName);
67 } 72 }
68 73
69 G4ENDFTapeRead::G4ENDFTapeRead(std::istringstr << 74 G4ENDFTapeRead::
70 G4FFGEnumeratio << 75 G4ENDFTapeRead( std::istringstream& dataStream,
71 G4FFGEnumeratio << 76 G4FFGEnumerations::YieldType WhichYield,
72 : /*Cause_(WhichCause),*/ << 77 G4FFGEnumerations::FissionCause /*WhichCause*/,
>> 78 G4int Verbosity )
>> 79 : /*Cause_(WhichCause),*/
73 Verbosity_(Verbosity), 80 Verbosity_(Verbosity),
74 YieldType_(WhichYield) 81 YieldType_(WhichYield)
75 { 82 {
76 // Initialize the class << 83 // Initialize the class
77 Initialize(dataStream); << 84 Initialize(dataStream);
78 } 85 }
79 86
80 void G4ENDFTapeRead::Initialize(const G4String << 87 void G4ENDFTapeRead::
>> 88 Initialize( G4String dataFile )
81 { 89 {
82 std::istringstream dataStream(std::ios::in); << 90 std::istringstream dataStream(std::ios::in);
83 G4ParticleHPManager::GetInstance()->GetDataS << 91 G4ParticleHPManager::GetInstance()->GetDataStream(dataFile, dataStream);
84 92
85 Initialize(dataStream); << 93 Initialize(dataStream);
86 } 94 }
87 95
88 void G4ENDFTapeRead::Initialize(std::istringst << 96 void G4ENDFTapeRead::
>> 97 Initialize( std::istringstream& dataStream )
>> 98 {
>> 99 G4FFG_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
>> 117 G4FFG_FUNCTIONLEAVE__
>> 118 }
>> 119
>> 120 G4double* G4ENDFTapeRead::
>> 121 G4GetEnergyGroupValues( void )
89 { 122 {
90 G4FFG_FUNCTIONENTER__ << 123 G4FFG_FUNCTIONENTER__
>> 124
>> 125 G4FFG_FUNCTIONLEAVE__
>> 126 return EnergyGroupValues_;
>> 127 }
91 128
92 EnergyGroups_ = 0; << 129 G4int G4ENDFTapeRead::
93 EnergyGroupValues_ = nullptr; << 130 G4GetNumberOfEnergyGroups( void )
>> 131 {
>> 132 G4FFG_FUNCTIONENTER__
94 133
95 YieldContainerTable_ = new G4TableTemplate<G << 134 G4FFG_FUNCTIONLEAVE__
>> 135 return EnergyGroups_;
>> 136 }
96 137
97 try { << 138 G4int G4ENDFTapeRead::
98 ReadInData(dataStream); << 139 G4GetNumberOfFissionProducts( void )
99 } << 140 {
100 catch (std::exception& e) { << 141 G4FFG_FUNCTIONENTER__
101 delete YieldContainerTable_; <<
102 142
103 G4FFG_FUNCTIONLEAVE__ << 143 G4int NumberOfElements = YieldContainerTable_->G4GetNumberOfElements();
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 144
447 dataIterator->second.first[currentEn << 145 G4FFG_FUNCTIONLEAVE__
448 dataIterator->second.second[currentE << 146 return NumberOfElements;
449 } << 147 }
450 else { << 148
451 // !!! Do not break since we need to << 149 G4ENDFYieldDataContainer* G4ENDFTapeRead::
452 // !!! entire data file for all poss << 150 G4GetYield( G4int WhichYield )
453 } << 151 {
454 } << 152 G4FFG_DATA_FUNCTIONENTER__
>> 153
>> 154 G4ENDFYieldDataContainer* Container = NULL;
>> 155 if(WhichYield >= 0 && WhichYield < YieldContainerTable_->G4GetNumberOfElements())
>> 156 {
>> 157 Container = YieldContainerTable_->G4GetContainer(WhichYield);
455 } 158 }
456 } <<
457 159
458 G4ENDFYieldDataContainer* NewDataContainer; << 160 G4FFG_DATA_FUNCTIONLEAVE__
459 EnergyGroups_ = (G4int)projectileEnergies.si << 161 return Container;
460 EnergyGroupValues_ = new G4double[EnergyGrou << 162 }
461 G4int NewProduct; << 163
462 G4FFGEnumerations::MetaState NewMetaState; << 164 void G4ENDFTapeRead::
463 auto NewYield = new G4double[EnergyGroups_]; << 165 G4SetVerbosity(G4int WhatVerbosity)
464 auto NewError = new G4double[EnergyGroups_]; << 166 {
465 << 167 G4FFG_FUNCTIONENTER__
466 for (G4int energyGroup = 0; energyGroup < En << 168
467 // Load the energy values << 169 this->Verbosity_ = WhatVerbosity;
468 EnergyGroupValues_[energyGroup] = projecti << 170
469 } << 171 G4FFG_FUNCTIONLEAVE__
470 << 172 }
471 // Load the data into the yield table << 173
472 for (dataIterator = intermediateData.begin() << 174 void G4ENDFTapeRead::
473 ++dataIterator) << 175 ReadInData( std::istringstream& dataStream )
474 { << 176 {
475 identifier = dataIterator->first; << 177 G4FFG_FUNCTIONENTER__
476 metastate = identifier % 10; << 178
477 switch (metastate) { << 179 // Check if the file exists
478 case 1: << 180 if(!dataStream.good())
479 NewMetaState = G4FFGEnumerations::META << 181 {
480 break; << 182 G4Exception("G4ENDFTapeRead::ReadInData()",
481 << 183 "Illegal file name",
482 case 2: << 184 JustWarning,
483 NewMetaState = G4FFGEnumerations::META << 185 "Fission product data not available");
484 break; << 186
485 << 187 // TODO create/use a specialized exception
486 default: << 188 G4FFG_FUNCTIONLEAVE__
487 G4Exception("G4ENDFTapeRead::ReadInDat << 189 throw std::exception();
488 "Unsupported metastable st <<
489 "the ground state"); <<
490 // Fall through <<
491 case 0: <<
492 NewMetaState = G4FFGEnumerations::GROU <<
493 break; <<
494 } 190 }
495 NewProduct = (identifier - metastate) / 10 <<
496 191
497 for (G4int energyGroup = 0; energyGroup < << 192 // Code to read in from a pure ENDF data tape.
498 if (energyGroup < (signed)dataIterator-> << 193 // Commented out while pre-formatted Geant4 ENDF data is being used
499 yield = dataIterator->second.first[ene << 194 // G4int CurrentEnergyGroup = -1;
500 error = dataIterator->second.second[en << 195 // std::vector< G4double > NewDoubleVector;
501 } << 196 // std::vector< G4double > EnergyPoints;
502 else { << 197 // std::vector< G4int > Product;
503 yield = 0.0; << 198 // std::vector< G4FFGEnumerations::MetaState > MetaState;
504 error = 0.0; << 199 // std::vector< std::vector< G4double > > Yield;
505 } << 200 // std::vector< std::vector< G4double > > Error;
>> 201 // G4String DataBlock;
>> 202 // 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 = 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 }
506 479
507 NewYield[energyGroup] = yield; << 480 G4ENDFYieldDataContainer* NewDataContainer;
508 NewError[energyGroup] = error; << 481 EnergyGroups_ = 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];
509 } 492 }
510 493
511 NewDataContainer = YieldContainerTable_->G << 494 // Load the data into the yield table
512 NewDataContainer->SetProduct(NewProduct); << 495 for(dataIterator = intermediateData.begin(); dataIterator != intermediateData.end(); ++dataIterator)
513 NewDataContainer->SetMetaState(NewMetaStat << 496 {
514 NewDataContainer->SetYieldProbability(NewY << 497 identifier = dataIterator->first;
515 NewDataContainer->SetYieldError(NewError); << 498 metastate = identifier % 10;
516 } << 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;
517 520
518 delete[] NewYield; << 521 for(G4int energyGroup = 0; energyGroup < EnergyGroups_; energyGroup++)
519 delete[] NewError; << 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 }
520 532
521 G4FFG_FUNCTIONLEAVE__ << 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
>> 547 G4FFG_FUNCTIONLEAVE__
522 } 548 }
523 549
524 G4ENDFTapeRead::~G4ENDFTapeRead() << 550 G4ENDFTapeRead::
>> 551 ~G4ENDFTapeRead( void )
525 { 552 {
526 G4FFG_FUNCTIONENTER__ << 553 G4FFG_FUNCTIONENTER__
527 554
528 delete[] EnergyGroupValues_; << 555 delete[] EnergyGroupValues_;
529 delete YieldContainerTable_; << 556 delete YieldContainerTable_;
530 557
531 G4FFG_FUNCTIONLEAVE__ << 558 G4FFG_FUNCTIONLEAVE__
532 } 559 }
>> 560
533 561