Start implementing 2D flux table for atmospheric fluxes

projects/distributions/private/primary/energy_direction/PrimaryEnergyDirectionDistribution.cxx

@@ -0,0 +1,31 @@
+#include "SIREN/distributions/primary/energy/PrimaryEnergyDirectionDistribution.h"
+
+#include <array>                                           // for array
+#include <string>                                          // for basic_string
+
+#include "SIREN/dataclasses/InteractionRecord.h"  // for Interactio...
+
+namespace siren {
+namespace distributions {
+
+//---------------
+// class PrimaryEnergyDirectionDistribution : PrimaryInjectionDistribution
+//---------------
+void PrimaryEnergyDirectionDistribution::Sample(
+        std::shared_ptr<siren::utilities::SIREN_random> rand,
+        std::shared_ptr<siren::detector::DetectorModel const> detector_model,
+        std::shared_ptr<siren::interactions::InteractionCollection const> interactions,
+        siren::dataclasses::PrimaryDistributionRecord & record) const {
+
+    std::pair<double,siren::math::Vector3D> energy_and_direction = SampleEnergyAndDirection(rand, detector_model, interactions, record);
+    record.SetEnergy(energy_and_direction.first);
+    record.SetDirection(energy_and_direction.second);
+}
+}
+
+std::vector<std::string> PrimaryEnergyDirectionDistribution::DensityVariables() const {
+    return std::vector<std::string>{"PrimaryEnergy","PrimaryDirection"};
+}
+
+} // namespace distributions
+} // namespace siren

projects/distributions/private/primary/energy_direction/Tabulated2DFluxDistribution.cxx

@@ -0,0 +1,271 @@
+#include "SIREN/distributions/primary/energy/Tabulated2DFluxDistribution.h"
+#include <array>                                           // for array
+#include <tuple>                                           // for tie, opera...
+#include <vector>                                          // for vector
+#include <fstream>                                         // for basic_istream
+#include <functional>                                      // for function
+
+#include "SIREN/dataclasses/InteractionRecord.h"  // for Interactio...
+#include "SIREN/distributions/Distributions.h"    // for InjectionD...
+#include "SIREN/utilities/Integration.h"          // for rombergInt...
+#include "SIREN/utilities/Interpolator.h"         // for TableData1D
+#include "SIREN/utilities/Random.h"               // for SIREN_random
+
+namespace siren { namespace interactions { class InteractionCollection; } }
+namespace siren { namespace detector { class DetectorModel; } }
+
+namespace siren {
+namespace distributions {
+namespace {
+    bool fexists(const std::string filename)
+    {
+            std::ifstream ifile(filename.c_str());
+            return (bool)ifile;
+    }
+}
+
+//---------------
+// class Tabulated2DFluxDistribution : PrimaryEnergyDirectionDistribution
+//---------------
+Tabulated2DFluxDistribution::Tabulated2DFluxDistribution() {}
+
+void Tabulated2DFluxDistribution::ComputeIntegral() {
+    std::function<double(double, double)> integrand = [&] (double x, double y) -> double {
+        return unnormed_pdf(x,y);
+    };
+    integral = siren::utilities::rombergIntegrate2D(integrand, energyMin, energyMax, zenithMin, zenithMax);
+}
+
+void Tabulated2DFluxDistribution::LoadFluxTable() {
+    if(fexists(fluxTableFilename)) {
+        std::ifstream in(fluxTableFilename.c_str());
+        std::string buf;
+        std::string::size_type pos;
+        siren::utilities::TableData2D<double> table_data;
+
+        while(std::getline(in, buf)) {
+            // Ignore comments and blank lines
+            if((pos = buf.find('#')) != std::string::npos)
+                buf.erase(pos);
+            const char* whitespace=" \n\r\t\v";
+            if((pos=buf.find_first_not_of(whitespace))!=0)
+                buf.erase(0,pos);
+            if(!buf.empty() && (pos=buf.find_last_not_of(whitespace))!=buf.size()-1)
+                buf.erase(pos+1);
+            if(buf.empty())
+                continue;
+
+            std::stringstream ss(buf);
+            double x, y, f;
+            ss >> x >> y >> f;
+            table_data.x.push_back(x);
+            table_data.y.push_back(y);
+            table_data.f.push_back(f);
+
+            energy_nodes.push_back(x);
+            zenith_nodes.push_back(y);
+        }
+        // If no physical are manually set, use first/last entry of table
+        if(not energy_bounds_set) {
+            energyMin = table_data.x[0];
+            energyMax = table_data.x[table_data.x.size()-1];
+        }
+        if(not zenith_bounds_set) {
+            zenithMin = table_data.y[0];
+            zenithMax = table_data.y[table_data.y.size()-1];
+        }
+        fluxTable = siren::utilities::Interpolator2D<double>(table_data);
+    } else {
+        throw std::runtime_error("Failed to open flux table file!");
+    }
+}
+
+void Tabulated2DFluxDistribution::LoadFluxTable(std::vector<double> & energies, std::vector<double> & zeniths, std::vector<double> & flux) {
+
+    assert(energies.size()==flux.size());
+    assert(zeniths.size()==flux.size());
+
+    siren::utilities::TableData2D<double> table_data;
+
+    table_data.x = energies;
+    table_data.y = zeniths;
+    table_data.f = flux;
+    energy_nodes = energies;
+    zenith_nodes = energies;
+
+    // If no physical are manually set, use first/last entry of table
+    if(not energy_bounds_set) {
+        energyMin = table_data.x[0];
+        energyMax = table_data.x[table_data.x.size()-1];
+    }
+    if(not zenith_bounds_set) {
+        zenithMin = table_data.y[0];
+        zenithMax = table_data.y[table_data.y.size()-1];
+    }
+    fluxTable = siren::utilities::Interpolator2D<double>(table_data);
+}
+
+double Tabulated2DFluxDistribution::unnormed_pdf(double energy, double zenith) const {
+    return fluxTable(energy, zenith);
+}
+
+double Tabulated2DFluxDistribution::SampleUnnormedPDF(double energy, double zenith) const {
+    return unnormed_pdf(energy, zenith);
+}
+
+double Tabulated2DFluxDistribution::GetIntegral() const {
+    return integral;
+}
+
+std::vector<double> Tabulated2DFluxDistribution::GetEnergyNodes() const {
+    return energy_nodes;
+}
+
+double Tabulated2DFluxDistribution::pdf(double energy, double zenith) const {
+    return unnormed_pdf(energy,zenith) / integral;
+}
+
+double Tabulated2DFluxDistribution::SamplePDF(double energy, double zenith) const {
+    return pdf(energy, zenith);
+}
+
+void Tabulated2DFluxDistribution::SetEnergyBounds(double eMin, double eMax) {
+    energyMin = eMin;
+    energyMax = eMax;
+    bounds_set = true;
+    ComputeIntegral();
+    ComputeCDF();
+}
+
+void Tabulated2DFluxDistribution::SetZenithBounds(double zMin, double zMax) {
+    zenithMin = zMin;
+    zenithMax = zMax;
+    bounds_set = true;
+    ComputeIntegral();
+    ComputeCDF();
+}
+
+Tabulated2DFluxDistribution::Tabulated2DFluxDistribution(std::string fluxTableFilename, bool has_physical_normalization)
+    : energy_bounds_set(false)
+    , zenith_bounds_set(false)
+    , fluxTableFilename(fluxTableFilename)
+{
+    LoadFluxTable();
+    ComputeIntegral();
+    if(has_physical_normalization)
+        SetNormalization(integral);
+}
+
+Tabulated2DFluxDistribution::Tabulated2DFluxDistribution(double energyMin, double energyMax, std::string fluxTableFilename, bool has_physical_normalization)
+    : energyMin(energyMin)
+    , energyMax(energyMax)
+    , energy_bounds_set(true)
+    , zenith_bounds_set(false)
+    , fluxTableFilename(fluxTableFilename)
+{
+    LoadFluxTable();
+    ComputeIntegral();
+    if(has_physical_normalization)
+        SetNormalization(integral);
+}
+
+Tabulated2DFluxDistribution::Tabulated2DFluxDistribution(double energyMin, double energyMax, double zenithMin, double zenithMax, std::string fluxTableFilename, bool has_physical_normalization)
+    : energyMin(energyMin)
+    , energyMax(energyMax)
+    , zenithMin(zenithMin)
+    , zenithMax(zenithMax)
+    , energy_bounds_set(true)
+    , zenith_bounds_set(true)
+    , fluxTableFilename(fluxTableFilename)
+{
+    LoadFluxTable();
+    ComputeIntegral();
+    if(has_physical_normalization)
+        SetNormalization(integral);
+}
+
+Tabulated2DFluxDistribution::Tabulated2DFluxDistribution(std::vector<double> energies, std::vector<double> zeniths, std::vector<double> flux, bool has_physical_normalization)
+    : energy_bounds_set(false)
+    , zenith_bounds_set(false)
+{
+    LoadFluxTable(energies,zeniths,flux);
+    ComputeIntegral();
+    if(has_physical_normalization)
+        SetNormalization(integral);
+}
+
+Tabulated2DFluxDistribution::Tabulated2DFluxDistribution(double energyMin, double energyMax, std::vector<double> energies, std::vector<double> zeniths, std::vector<double> flux, bool has_physical_normalization)
+    : energyMin(energyMin)
+    , energyMax(energyMax)
+    , energy_bounds_set(true)
+    , zenith_bounds_set(false)
+{
+    LoadFluxTable(energies,zeniths,flux);
+    ComputeIntegral();
+    if(has_physical_normalization)
+        SetNormalization(integral);
+}
+
+Tabulated2DFluxDistribution::Tabulated2DFluxDistribution(double energyMin, double energyMax, double zenithMin, double zenithMax, std::vector<double> energies, std::vector<double> zeniths, std::vector<double> flux, bool has_physical_normalization)
+    : energyMin(energyMin)
+    , energyMax(energyMax)
+    : zenithMin(zenithMin)
+    , zenithMax(zenithMax)
+    , energy_bounds_set(true)
+    , zenith_bounds_set(true)
+{
+    LoadFluxTable(energies,zeniths,flux);
+    ComputeIntegral();
+    if(has_physical_normalization)
+        SetNormalization(integral);
+}
+
+double Tabulated2DFluxDistribution::SampleEnergy(std::shared_ptr<siren::utilities::SIREN_random> rand, std::shared_ptr<siren::detector::DetectorModel const> detector_model, std::shared_ptr<siren::interactions::InteractionCollection const> interactions, siren::dataclasses::PrimaryDistributionRecord & record) const {
+    // TODO: Use metropolis hastings
+    return 0;
+}
+
+
+double Tabulated2DFluxDistribution::GenerationProbability(std::shared_ptr<siren::detector::DetectorModel const> detector_model, std::shared_ptr<siren::interactions::InteractionCollection const> interactions, siren::dataclasses::InteractionRecord const & record) const {
+    double const & energy = record.primary_momentum[0];
+    if(energy < energyMin or energy > energyMax)
+        return 0.0;
+    // TODO: check zenith
+    else if(energy < energyMin or energy > energyMax)
+        return 0.0;
+    else
+        return pdf(energy);
+}
+
+std::string Tabulated2DFluxDistribution::Name() const {
+    return "Tabulated2DFluxDistribution";
+}
+
+std::shared_ptr<PrimaryInjectionDistribution> Tabulated2DFluxDistribution::clone() const {
+    return std::shared_ptr<PrimaryInjectionDistribution>(new Tabulated2DFluxDistribution(*this));
+}
+
+bool Tabulated2DFluxDistribution::equal(WeightableDistribution const & other) const {
+    const Tabulated2DFluxDistribution* x = dynamic_cast<const Tabulated2DFluxDistribution*>(&other);
+
+    if(!x)
+        return false;
+    else
+        return
+            std::tie(energyMin, energyMax, zenithMin, zenithMax, fluxTable)
+            ==
+            std::tie(x->energyMin, x->energyMax, x->zenithMin, x->zenithMax, x->fluxTable);
+}
+
+bool Tabulated2DFluxDistribution::less(WeightableDistribution const & other) const {
+    const Tabulated2DFluxDistribution* x = dynamic_cast<const Tabulated2DFluxDistribution*>(&other);
+    return
+        std::tie(energyMin, energyMax, zenithMin, zenithMax, fluxTable)
+        <
+        std::tie(x->energyMin, x->energyMax, x->zenithMin, x->zenithMax, x->fluxTable);
+}
+
+} // namespace distributions
+
+} // namespace siren
+

projects/distributions/public/SIREN/distributions/primary/energy_direction/PrimaryEnergyDirectionDistribution.h

@@ -0,0 +1,71 @@
+#pragma once
+#ifndef SIREN_PrimaryEnergyDirectionDistribution_H
+#define SIREN_PrimaryEnergyDirectionDistribution_H
+
+#include <memory>                                        // for shared_ptr
+#include <string>                                        // for string
+#include <vector>                                        // for vector
+#include <cstdint>                                       // for uint32_t
+#include <stdexcept>                                     // for runtime_error
+
+#include <cereal/access.hpp>
+#include <cereal/types/polymorphic.hpp>
+#include <cereal/types/base_class.hpp>
+#include <cereal/types/utility.hpp>
+
+#include "SIREN/distributions/Distributions.h"  // for WeightableDi...
+#include "SIREN/math/Vector3D.h"  // for Vector3D
+
+namespace siren { namespace interactions { class InteractionCollection; } }
+namespace siren { namespace dataclasses { class InteractionRecord; } }
+namespace siren { namespace detector { class DetectorModel; } }
+namespace siren { namespace utilities { class SIREN_random; } }
+namespace cereal { class access; }
+
+namespace siren {
+namespace distributions {
+
+class PrimaryEnergyDirectionDistribution : virtual public PrimaryInjectionDistribution, virtual public PhysicallyNormalizedDistribution {
+friend cereal::access;
+public:
+    virtual ~PrimaryEnergyDirectionDistribution() {};
+private:
+    virtual std::pair<double,siren::math::Vector3D> SampleEnergyAndDirection(std::shared_ptr<siren::utilities::SIREN_random> rand, std::shared_ptr<siren::detector::DetectorModel const> detector_model, std::shared_ptr<siren::interactions::InteractionCollection const> interactions, siren::dataclasses::PrimaryDistributionRecord & record) const = 0;
+public:
+    void Sample(std::shared_ptr<siren::utilities::SIREN_random> rand, std::shared_ptr<siren::detector::DetectorModel const> detector_model, std::shared_ptr<siren::interactions::InteractionCollection const> interactions, siren::dataclasses::PrimaryDistributionRecord & record) const override;
+    virtual double GenerationProbability(std::shared_ptr<siren::detector::DetectorModel const> detector_model, std::shared_ptr<siren::interactions::InteractionCollection const> interactions, siren::dataclasses::InteractionRecord const & record) const override = 0;
+    virtual std::vector<std::string> DensityVariables() const override;
+    virtual std::string Name() const override = 0;
+    virtual std::shared_ptr<PrimaryInjectionDistribution> clone() const override = 0;
+    template<typename Archive>
+    void save(Archive & archive, std::uint32_t const version) const {
+        if(version == 0) {
+            archive(cereal::virtual_base_class<PrimaryInjectionDistribution>(this));
+            archive(cereal::virtual_base_class<PhysicallyNormalizedDistribution>(this));
+        } else {
+            throw std::runtime_error("PrimaryEnergyDirectionDistribution only supports version <= 0!");
+        }
+    }
+    template<typename Archive>
+    void load(Archive & archive, std::uint32_t const version) {
+        if(version == 0) {
+            archive(cereal::virtual_base_class<PrimaryInjectionDistribution>(this));
+            archive(cereal::virtual_base_class<PhysicallyNormalizedDistribution>(this));
+        } else {
+            throw std::runtime_error("PrimaryEnergyDirectionDistribution only supports version <= 0!");
+        }
+    }
+protected:
+    virtual bool equal(WeightableDistribution const & distribution) const override = 0;
+    virtual bool less(WeightableDistribution const & distribution) const override = 0;
+};
+
+} // namespace distributions
+} // namespace siren
+
+CEREAL_CLASS_VERSION(siren::distributions::PrimaryEnergyDirectionDistribution, 0);
+CEREAL_REGISTER_TYPE(siren::distributions::PrimaryEnergyDirectionDistribution);
+CEREAL_REGISTER_POLYMORPHIC_RELATION(siren::distributions::PrimaryInjectionDistribution, siren::distributions::PrimaryEnergyDirectionDistribution);
+CEREAL_REGISTER_POLYMORPHIC_RELATION(siren::distributions::PhysicallyNormalizedDistribution, siren::distributions::PrimaryEnergyDirectionDistribution);
+
+#endif // SIREN_PrimaryEnergyDirectionDistribution_H