Commit e1cafb65 authored by LEFEBVREJP email's avatar LEFEBVREJP email
Browse files

Merge branch 'f71togapi' into 'master'

F71 to GAPI

See merge request !120
parents debfc7cb 56778b2f
Pipeline #159242 passed with stages
in 18 minutes and 36 seconds
......@@ -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)
......
#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
#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_ */
#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;
}
......@@ -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)
......
#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);
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment