emitdb.cc 6.27 KB
Newer Older
1
2
3
4
5
#include <fstream>

#include "radixcore/stringfunctions.hh"
#include "radixcore/system.hh"

6
7
8
9
10
11
#include "radixio/emitdb.hh"

namespace radix
{
EmitDb::EmitDb()
{
12
  // load_mpdkxgam(EmitDb::XGAMDATA4);
13
  // load_mpdkxgam(EmitDb::XGAMDATA3);
14
  // load_ensdfbeta(EmitDb::BETADATA);
15
16
17
18
19
20
}

EmitDb::EmitDb(const EmitDb& orig)
{
  mIsotope      = orig.mIsotope;
  mIsotopeTitle = orig.mIsotopeTitle;
21
22
}

23
void XGamDb::load(const std::string& file, EmitDb& database)
24
25
{
  std::string revStr = radix::split_string(".", file)[1];
26
  int revNum         = std::atoi(revStr.substr(revStr.size() - 1, 1).c_str());
27
  int indexThreshold = 70, indexIterator = 20;
28
29
30
31
  int linecount = 0;
  if (revNum == 3)
  {
    indexThreshold = 70;
32
    indexIterator  = 20;
33
34
35
36
  }
  else if (revNum == 4)
  {
    indexThreshold = 77;
37
    indexIterator  = 22;
38
39
40
41
  }
  std::ifstream xgamFile(file);
  std::string line, prev;
  std::vector<float> energies, intensities;
42
43
  int zaid          = -1;
  int totalLines    = 0;
44
  float totalEnergy = 0.f;
45
  auto& isotopes    = database.isotopes();
46
47
48
49
50
51
52
53
54
  while (std::getline(xgamFile, line))
  {
    linecount += 1;
    int thisZaid = std::atoi(line.substr(1, 6).c_str());
    if (thisZaid != 0)
    {
      if (zaid != -1 && energies.size() > 0)
      {
        EmitDBEntry emitEntry;
55
56
57
58
59
60
61
        emitEntry.z           = int(zaid / 10000);
        emitEntry.a           = int((zaid % 10000) / 10);
        emitEntry.m           = int(zaid % 10);
        emitEntry.energy      = totalEnergy;
        emitEntry.energies    = energies;
        emitEntry.intensities = intensities;
        isotopes.insert(std::make_pair(zaid, emitEntry));
62
63
64
        energies.clear();
        intensities.clear();
      }
65
66
      zaid = thisZaid;
      // totalLines  = std::atoi(line.substr(9, 12).c_str());
67
68
69
70
      totalEnergy = std::atof(line.substr(60, 10).c_str());
    }
    else
    {
71
72
      int index      = 3;
      int jumpLength = 10;
73
74
      if (revNum == 3)
      {
75
        index      = 10;
76
77
78
79
        jumpLength = 10;
      }
      else
      {
80
        index      = 10;
81
82
83
84
85
        jumpLength = 11;
      }
      while (index < indexThreshold && index < line.size())
      {
        float energy = std::atof(line.substr(index, 10).c_str());
86
87
        float intensity =
            std::atof(line.substr(index + jumpLength, 10).c_str());
88
89
90
91
92
93
94
95
        energies.push_back(energy);
        intensities.push_back(intensity);
        index += indexIterator;
      }
    }
    prev = line;
  }
  EmitDBEntry emitEntry;
96
97
98
99
100
101
102
  emitEntry.z           = int(zaid / 10000);
  emitEntry.a           = int((zaid % 10000) / 10);
  emitEntry.m           = int(zaid % 10);
  emitEntry.energy      = totalEnergy;
  emitEntry.energies    = energies;
  emitEntry.intensities = intensities;
  isotopes.insert(std::make_pair(zaid, emitEntry));
103
104
  energies.clear();
  intensities.clear();
105
}
106

107
void EnsdfBetaDb::load(const std::string& file, EmitDb& database)
108
109
110
111
{
  std::ifstream betaFile(file);
  std::string line, prev;
  int zaidCount = 0;
112
  int zaid      = -1;
113
114
  int linecount = 0;
  std::vector<float> energies, intensities;
115
  auto& isotopes = database.isotopes();
116
117
118
119
120
121
122
  while (std::getline(betaFile, line))
  {
    linecount += 1;
    if (line.substr(0, 2) == "  ")
    {
      if (line.size() > 10)
      {
123
124
        int zaidl = std::atoi(line.substr(2, 6).c_str());
        if (zaid != zaidl)
125
        {
126
          zaid = zaidl;
127
        }
128
        zaidCount++;
129
        float intensity = std::atof(line.substr(10, 8).c_str());
130
        float energy    = std::atof(line.substr(19, 9).c_str());
131
132
133
134
135
136
137
138
139
140
141
142
143
        energies.push_back(energy);
        intensities.push_back(intensity);
      }
      else
      {
        if (energies.size() > 0 && zaid != -1)
        {
          float totalEnergy = 0.0f;
          for (size_t i = 0; i < intensities.size(); ++i)
          {
            totalEnergy += (intensities[i] * energies[i]);
          }
          EmitDBEntry emitEntry;
144
145
146
147
148
149
          emitEntry.z           = int(zaid / 10000);
          emitEntry.a           = int((zaid % 10000) / 10);
          emitEntry.m           = int(zaid % 10);
          emitEntry.energy      = totalEnergy;
          emitEntry.energies    = energies;
          emitEntry.intensities = intensities;
150
          const auto& isOutside =
151
              isotopes.insert(std::make_pair(zaid, emitEntry));
152
153
          if (!isOutside.second)
          {
154
155
156
            isOutside.first->second.energy      = totalEnergy;
            isOutside.first->second.energies    = energies;
            isOutside.first->second.intensities = intensities;
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
          }
          zaidCount = 0;
          energies.clear();
          intensities.clear();
        }
      }
    }
    prev = line;
  }
  if (energies.size() > 0)
  {
    float totalEnergy = 0.0f;
    for (size_t i = 0; i < intensities.size(); ++i)
    {
      totalEnergy += (intensities[i] * energies[i]);
    }
    EmitDBEntry emitEntry;
174
175
176
177
178
179
180
    emitEntry.z           = int(zaid / 10000);
    emitEntry.a           = int((zaid % 10000) / 10);
    emitEntry.m           = int(zaid % 10);
    emitEntry.energy      = totalEnergy;
    emitEntry.energies    = energies;
    emitEntry.intensities = intensities;
    const auto& isOutside = isotopes.insert(std::make_pair(zaid, emitEntry));
181
182
    if (!isOutside.second)
    {
183
184
185
      isOutside.first->second.energy      = totalEnergy;
      isOutside.first->second.energies    = energies;
      isOutside.first->second.intensities = intensities;
186
187
188
189
190
191
    }
    energies.clear();
    intensities.clear();
  }
}

192
193
194
195
196
197
198
199
200
201
const std::string& EmitDb::isotope_title() const { return mIsotopeTitle; }
void EmitDb::set_isotope_title(const std::string& title)
{
  mIsotopeTitle = title;
}
const std::unordered_map<int, EmitDBEntry>& EmitDb::isotopes() const
{
  return mIsotope;
}
std::unordered_map<int, EmitDBEntry>& EmitDb::isotopes() { return mIsotope; }
202
const EmitDBEntry& EmitDb::emitter(int zaid) const
203
204
205
206
{
  const auto& it = mIsotope.find(zaid);
  return it->second;
}
207
std::vector<std::vector<float>> EmitDb::energy_intensities(const int zaid) const
208
209
{
  const auto& it                         = mIsotope.find(zaid);
210
211
  const std::vector<float>& energies     = it->second.energies;
  const std::vector<float>& intensities  = it->second.intensities;
212
213
214
215
  std::vector<std::vector<float>> paired = {energies, intensities};
  return paired;
}

216
217
218
219
220
221
void EmitDb::operator=(const EmitDb& orig)
{
  mIsotope      = orig.mIsotope;
  mIsotopeTitle = orig.mIsotopeTitle;
}

222
}  // namespace radix