F71 to GAPI utility

radixio/CMakeLists.txt

@@ -11,6 +11,7 @@ SET(SOURCE
    decaydb.cc
    gfsfile.cc
    f71stream.cc
+   f71togapi.cc
    hcdump.cc
    hysplitcdump.cc
    spectrum.cc
@@ -28,6 +29,7 @@ SET(HEADERS
    gfsfile.hh
    f71stream.hh
    f71stream.i.hh
+   f71togapi.hh
    hcdump.hh
    hysplitcdump.hh
    hysplitcdump.i.hh
@@ -43,6 +45,15 @@ TRIBITS_ADD_LIBRARY(radixio
   SOURCES ${SOURCE}
   NOINSTALLHEADERS ${HEADERS}
 )
+
+#
+# Add executables
+TRIBITS_ADD_EXECUTABLE(radixf71togapi
+    NOEXEPREFIX NOEXESUFFIX
+    SOURCES f71togapi_main.cc
+    INSTALLABLE
+    )
+
 #
 # Add testing directory
 TRIBITS_ADD_TEST_DIRECTORIES(tests)

radixio/f71togapi.cc

+#include "radixio/f71togapi.hh"
+
+#include "radixio/decaydb.hh"
+#include "radixio/zaid.hh"
+
+namespace radix
+{
+void F71toGapi::to_gapi_json(const F71Case& fcase, float detector_height,
+                             float detector_distance, float live_time,
+                             DataObject& obj, float threshold)
+{
+  //
+  //
+  // calculate activity (curies) for each nuclide
+  const std::vector<int>& light_nucs = fcase.lightElementNuclides();
+  const std::vector<int>& ac_nucs    = fcase.actinideNuclides();
+  const std::vector<int>& fp_nucs    = fcase.fissionProductNuclides();
+
+  const std::vector<float>& light_abs = fcase.lightElementNuclideAbundances();
+  const std::vector<float>& ac_abs    = fcase.actinideNuclideAbundances();
+  const std::vector<float>& fp_abs    = fcase.fissionProductNuclideAbundances();
+  // collapse nuclides to a single set with all three abundances
+  std::unordered_map<int, std::tuple<float, float, float>> master;
+  // process light nuclides
+  for (size_t i = 0; i < light_nucs.size(); ++i)
+  {
+    int zaid  = light_nucs[i];
+    float abs = light_abs[i];
+    auto it   = master.find(zaid);
+    if (it == master.end())
+    {
+      master.insert(std::make_pair(zaid, std::make_tuple(abs, 0.f, 0.f)));
+    }
+    // we won't have this nuclide in the master list before
+  }
+  // process actinide nuclides
+  for (size_t i = 0; i < ac_nucs.size(); ++i)
+  {
+    int zaid  = ac_nucs[i];
+    float abs = ac_abs[i];
+    auto it   = master.find(zaid);
+    if (it == master.end())
+    {
+      master.insert(std::make_pair(zaid, std::make_tuple(0.f, abs, 0.f)));
+    }
+    else
+    {
+      std::get<1>(it->second) = abs;
+    }
+  }
+  // process fission product nuclides
+  for (size_t i = 0; i < fp_nucs.size(); ++i)
+  {
+    int zaid  = fp_nucs[i];
+    float abs = fp_abs[i];
+    auto it   = master.find(zaid);
+    if (it == master.end())
+    {
+      master.insert(std::make_pair(zaid, std::make_tuple(0.f, 0.f, abs)));
+    }
+    else
+    {
+      std::get<2>(it->second) = abs;
+    }
+  }
+  DecayDb decaydb;
+  std::vector<std::pair<int, float>> isotope_activity;
+  for (const auto& it : master)
+  {
+    int zaid        = it.first;
+    float light_abs = std::get<0>(it.second);
+    float ac_abs    = std::get<1>(it.second);
+    float fp_abs    = std::get<2>(it.second);
+    float activity =
+        decaydb.calculate_activity(zaid, light_abs, ac_abs, fp_abs);
+    radix_line("Calculated activity(" << activity << ") for " << zaid << "["
+                                      << light_abs << ", " << ac_abs << ", "
+                                      << fp_abs << "]");
+    if (threshold < activity) isotope_activity.push_back({zaid, activity});
+  }
+  //
+  // Begin DataObject creation
+  obj["detector height"]   = detector_height;
+  obj["detector distance"] = detector_distance;
+  obj["live time"]         = live_time;
+  DataArray sources;
+  for (const auto& isotope : isotope_activity)
+  {
+    DataObject source;
+    source["name"]     = zaid_symbol(isotope.first);
+    source["activity"] = isotope.second;
+    source["age"]      = 0;
+    DataArray shields;  // no shields for now
+    source["shields"] = shields;
+    sources.push_back(source);
+  }
+  obj["sources"] = sources;
+}
+}  // namespace radix

radixio/f71togapi.hh

+#ifndef RADIX_RADIXIO_F71TOGAPI_HH_
+#define RADIX_RADIXIO_F71TOGAPI_HH_
+
+#include "radixcore/value.hh"
+#include "radixcore/visibility.hh"
+
+#include "radixio/f71stream.hh"
+
+namespace radix
+{
+class RADIX_PUBLIC F71toGapi
+{
+ public:
+  /**
+   * @brief to_gapi_json Translate F71 isotopics to gadras api isotope and
+   * activity json object
+   * @param fcase F71Case
+   * @param detector_height height of the detector in cm above the ground
+   * @param detector_distance distance in cm to source
+   * @param live_time measurement time in seconds
+   * @param obj DataObject to construct
+   * @param threshold in curies
+   */
+  static void to_gapi_json(const F71Case& fcase, float detector_height,
+                           float detector_distance, float live_time,
+                           DataObject& obj, float threshold = 1e-6);
+};
+}  // namespace radix
+
+#endif /* RADIX_RADIXIO_F71TOGAPI_HH_ */

radixio/f71togapi_main.cc

+#include "radixcore/commandline.hh"
+#include "radixcore/system.hh"
+#include "radixio/f71togapi.hh"
+
+#include <string>
+#include <vector>
+
+using namespace radix;
+int main(int argc, char** argv)
+{
+  CommandLine commandline(argc, argv);
+  commandline.setDescription("RADIX F71 file to GADRAS API (GAPI) JSON input.");
+  commandline.addOption("f71", "SCALE6.1 F71 formatted binary file.", true);
+  commandline.addOption("o", "JSON output file. Defaults to f71 basename.json");
+  commandline.addOption("h", "Detector height (cm)", true);
+  commandline.addOption("d", "Detector distance (cm)", true);
+  commandline.addOption("l", "Live time (seconds)", true);
+  commandline.addOption("i", "F71 case index. Defaults to 0");
+
+  {
+    std::vector<std::string> errors;
+    if (!commandline.validate(errors))
+    {
+      for (const auto& e : errors)
+      {
+        std::cerr << e << std::endl;
+      }
+      commandline.help(std::cerr);
+      return 1;
+    }
+  }
+  std::string f71file     = commandline.get<std::string>("f71");
+  std::string output_file = commandline.get<std::string>("o", "");
+  float detector_height   = commandline.get<float>("h");
+  float detector_distance = commandline.get<float>("d");
+  float live_time         = commandline.get<float>("l");
+  int index               = commandline.get<int>("i", 0);
+  if (output_file.empty())
+  {
+    output_file = to_native_path(dirname(f71file) + "/" +
+                                 basename(f71file, true) + ".json");
+    std::cout << "**Notice: No output file specified. Defaulting to '"
+              << output_file << "'" << std::endl;
+  }
+  if (live_time < 1)
+    std::cerr << "***Error: Live time can not be less than 1." << std::endl;
+  DataObject obj;
+  F71Stream stream(f71file);
+  bool result   = true;
+  size_t findex = 0;
+  bool success  = false;
+  while (result)
+  {
+    F71Case single_case;
+    try
+    {
+      result = stream.read_case(single_case);
+    }
+    catch (const std::exception& e)
+    {
+      std::cerr << "***Error: " << e.what() << std::endl;
+      break;
+    }
+
+    if (result && findex == index)
+    {
+      F71toGapi::to_gapi_json(single_case, detector_height, detector_distance,
+                              live_time, obj);
+      std::ofstream ofs(output_file);
+      obj.pack_json(ofs);
+      ofs.close();
+      success = true;
+      break;
+    }
+    else
+    {
+      std::cerr << "***Error: Failed to read case " << findex << std::endl;
+      break;
+    }
+    findex++;
+  }
+  return (success) ? 0 : 1;
+}

radixio/tests/CMakeLists.txt

@@ -4,6 +4,7 @@ ADD_GOOGLE_TEST(tstCSVFile.cc NP 1)
 ADD_GOOGLE_TEST(tstCFGFile.cc NP 1)
 ADD_GOOGLE_TEST(tstEafstream.cc NP 1)
 ADD_GOOGLE_TEST(tstF71Stream.cc NP 1)
+ADD_GOOGLE_TEST(tstF71toGapi.cc NP 1)
 ADD_GOOGLE_TEST(tstDecayDb.cc NP 1)
 ADD_GOOGLE_TEST(tstHysplitCDump.cc NP 1)
 ADD_GOOGLE_TEST(tstSpectrum.cc NP 1)

radixio/tests/tstF71toGapi.cc

+
+#include "gtest/gtest.h"
+#include "radixcore/system.hh"
+#include "radixio/f71togapi.hh"
+
+#include <fstream>
+using namespace radix;
+TEST(Radixio, F71toGapi)
+{
+  // test reading a scale 6.1 formated f71 file from DELFIC
+  F71Stream stream(radix::to_native_path(
+      std::string(dirname(__FILE__) + "/data/test.s61.f71")));
+  std::vector<F71Case> cases;
+  bool result = true;
+  while (result)
+  {
+    F71Case single_case;
+    result = stream.read_case(single_case);
+    if (result) cases.emplace_back(single_case);
+  }
+  ASSERT_TRUE(cases.size() == 1);
+  F71Case& first = cases[0];
+  DataObject obj;
+  F71toGapi::to_gapi_json(first, 100.f, 1000.f, 1, obj);
+  EXPECT_EQ(100., obj["detector height"].to_double());
+  EXPECT_EQ(1000., obj["detector distance"].to_double());
+  EXPECT_EQ(1., obj["live time"].to_double());
+  ASSERT_TRUE(obj.contains("sources"));
+  const DataArray& sources = obj["sources"].as_array();
+  EXPECT_EQ(231, sources.size());
+  float total_activity = 0.;
+  float blessed_total  = 253.6866455078125;
+  for (size_t i = 0; i < sources.size(); ++i)
+  {
+    const DataObject& source = sources[i].as_object();
+    std::string name         = source["name"].to_string();
+    float activity           = float(source["activity"].to_double());
+    std::cout << name << " " << activity << " curies" << std::endl;
+    total_activity += activity;
+  }
+  EXPECT_NEAR(blessed_total, total_activity, 1e-3);
+}