source_sampling.h

Go to the documentation of this file.

#ifndef PYNE_6OR6BJURKJHHTOFWXO2VMQM5EY
#define PYNE_6OR6BJURKJHHTOFWXO2VMQM5EY

#include <assert.h>
#include <iostream>
#include <fstream>
#include <stdio.h>
#include <stdlib.h>
#include <vector>
#include <stdexcept>
#include <sstream>
#include <string>
#include <map>

#include "moab/Range.hpp"
#include "moab/Core.hpp"
#ifndef PYNE_IS_AMALGAMATED
#include "measure.h"
#endif
#include "moab/CartVect.hpp"

#ifdef __cplusplus
extern "C" {
#endif

namespace pyne {

  void sampling_setup_(int* mode, int* cell_list_size);
  void particle_birth_(double* rands,
                            double* x,
                            double* y,
                            double* z,
                            double* e,
                            double* w,
                            int* cell_list);
  std::vector<double> read_e_bounds(std::string e_bounds_file);

  struct edge_points {
    moab::CartVect o_point;
    moab::CartVect x_vec;
    moab::CartVect y_vec;
    moab::CartVect z_vec;
  };

  class AliasTable {
  public:
    AliasTable(std::vector<double> p);
    int sample_pdf(double rand1, double rand2);
    ~AliasTable(){};
    int n; 
    std::vector<double> prob; 
    std::vector<int> alias; 
  };

  // class Source particle
  class SourceParticle {
    public:
    SourceParticle();
    SourceParticle(double x,
                   double y,
                   double z,
                   double e,
                   double w,
                   std::vector<int> cell_list);
    ~SourceParticle();

    double get_x() {return x;};
    double get_y() {return y;};
    double get_z() {return z;};
    double get_e() {return e;};
    double get_w() {return w;};
    std::vector<int> get_cell_list() {return cell_list;};
    private:
    double x; // x coordinate
    double y; // y coordinate
    double z; // z coordinate
    double e; // energy
    double w; // weight
    std::vector<int> cell_list;
  };

  enum BiasMode {USER, ANALOG, UNIFORM};
  enum MeshMode {VOXEL, SUBVOXEL, TET};

  class Sampler {
  public:
    Sampler(std::string filename,
            std::string src_tag_name,
            std::vector<double> e_bounds,
            bool uniform);
    Sampler(std::string filename,
            std::string src_tag_name,
            std::vector<double> e_bounds,
            std::string bias_tag_name);
    Sampler(std::string filename,
            std::map<std::string, std::string> tag_names,
            std::vector<double> e_bounds,
            int mode);

    pyne::SourceParticle particle_birth(std::vector<double> rands);

    int get_cell_list_size();

    ~Sampler() {
      delete mesh;
      delete at;
    };

  // member variables
  private:
    // problem parameters
    std::string filename; 
    std::string src_tag_name; 
    std::string bias_tag_name; 
    std::string cell_number_tag_name; 
    std::string cell_fracs_tag_name; 
    std::map<std::string, std::string> tag_names; 
    std::vector<double> e_bounds;  
    int num_e_groups; 
    int num_bias_groups; 
    int max_num_cells; 
    // Equivalent cell number in photon source.
    // For voxel R2S, p_src_num_cells = 1
    // For sub-voxel R2S, p_src_num_cells = max_num_cells
    // For unstructured R2S, p_src_num_cells = 1.
    int p_src_num_cells;
    int cell_list_size;
    bool has_cell_fracs;
    BiasMode bias_mode; 
    MeshMode mesh_mode; 
    int mode; 
    // mesh
    moab::Interface* mesh; 
    int num_ves; 
    moab::EntityType ve_type; 
    int verts_per_ve; 
    // sampling
    std::vector<edge_points> all_edge_points; 
    std::vector<double> biased_weights; 
    std::vector<int> cell_number; 
    std::vector<int> cell_list; 
    std::vector<double> cell_fracs; 
    AliasTable* at; 

  // member functions
  private:
    // instantiation
    void setup();
    void mesh_geom_data(moab::Range ves, std::vector<double> &volumes);
    void mesh_tag_data(moab::Range ves, const std::vector<double> volumes);
    // select birth parameters
    moab::CartVect sample_xyz(int ve_idx, std::vector<double> rands);
    double sample_e(int e_idx, double rand);
    double sample_w(int pdf_idx);
    // helper functions
    void normalize_pdf(std::vector<double> & pdf);
    int num_groups(moab::Tag tag);
    std::vector<double> read_bias_pdf(moab::Range ves, std::vector<double> volumes,
                                      std::vector<double> pdf);
    // Get max_num_cells
    int get_max_num_cells(moab::Tag cell_fracs_tag);
    // get has_cell_fracs
    bool check_cell_fracs(moab::Tag cell_fracs_tag);
  };
} //end namespace pyne

#ifdef __cplusplus
} // extern "C"
#endif

#endif // PYNE_6OR6BJURKJHHTOFWXO2VMQM5EY