ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ENDFTapeRead.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ENDFTapeRead.cc
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 
47 G4ENDFTapeRead::
48 G4ENDFTapeRead( 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 
60 G4ENDFTapeRead::
61 G4ENDFTapeRead( 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 
74 G4ENDFTapeRead::
75 G4ENDFTapeRead( 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 
87 void G4ENDFTapeRead::
88 Initialize( G4String dataFile )
89 {
90  std::istringstream dataStream(std::ios::in);
91  G4ParticleHPManager::GetInstance()->GetDataStream(dataFile, dataStream);
92 
93  Initialize(dataStream);
94 }
95 
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 )
122 {
123 G4FFG_FUNCTIONENTER__
124 
125 G4FFG_FUNCTIONLEAVE__
126  return EnergyGroupValues_;
127 }
128 
129 G4int G4ENDFTapeRead::
130 G4GetNumberOfEnergyGroups( void )
131 {
132 G4FFG_FUNCTIONENTER__
133 
134 G4FFG_FUNCTIONLEAVE__
135  return EnergyGroups_;
136 }
137 
138 G4int G4ENDFTapeRead::
139 G4GetNumberOfFissionProducts( void )
140 {
141 G4FFG_FUNCTIONENTER__
142 
143  G4int NumberOfElements = YieldContainerTable_->G4GetNumberOfElements();
144 
145 G4FFG_FUNCTIONLEAVE__
146  return NumberOfElements;
147 }
148 
149 G4ENDFYieldDataContainer* G4ENDFTapeRead::
150 G4GetYield( G4int WhichYield )
151 {
152 G4FFG_DATA_FUNCTIONENTER__
153 
154  G4ENDFYieldDataContainer* Container = NULL;
155  if(WhichYield >= 0 && WhichYield < YieldContainerTable_->G4GetNumberOfElements())
156  {
157  Container = YieldContainerTable_->G4GetContainer(WhichYield);
158  }
159 
160 G4FFG_DATA_FUNCTIONLEAVE__
161  return Container;
162 }
163 
164 void G4ENDFTapeRead::
165 G4SetVerbosity(G4int WhatVerbosity)
166 {
167 G4FFG_FUNCTIONENTER__
168 
169  this->Verbosity_ = WhatVerbosity;
170 
171 G4FFG_FUNCTIONLEAVE__
172 }
173 
174 void G4ENDFTapeRead::
175 ReadInData( std::istringstream& dataStream )
176 {
177 G4FFG_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 // 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  }
479 
480  G4ENDFYieldDataContainer* NewDataContainer;
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];
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 
547 G4FFG_FUNCTIONLEAVE__
548 }
549 
550 G4ENDFTapeRead::
551 ~G4ENDFTapeRead( void )
552 {
553 G4FFG_FUNCTIONENTER__
554 
555  delete[] EnergyGroupValues_;
556  delete YieldContainerTable_;
557 
558 G4FFG_FUNCTIONLEAVE__
559 }
560