ResponseTable.hh

//
// This file is part of MARLEY (Model of Argon Reaction Low Energy Yields)
//
// MARLEY is free software: you can redistribute it and/or modify it under the
// terms of version 3 of the GNU General Public License as published by the
// Free Software Foundation.
//
// For the full text of the license please see COPYING or
// visit http://opensource.org/licenses/GPL-3.0
//
// Please respect the MCnet academic usage guidelines. See GUIDELINES
// or visit https://www.montecarlonet.org/GUIDELINES for details.
 
#pragma once
 
// Standard library includes
#include <fstream>
#include <functional>
#include <memory>
#include <vector>
 
// MARLEY includes
#include "marley/NuclearResponses.hh"
 
namespace marley {

  class ResponseTable {
    public:
 
      ResponseTable() {}
 
      ResponseTable( std::shared_ptr< std::vector<double> >& wvec,
        std::shared_ptr< std::vector<double> >& qvec,
        std::shared_ptr< std::vector<NuclearResponses> >& resp )
        : wvec_( wvec ), qvec_( qvec ), responses_( resp ) {}

      inline double w_min() const { return wvec_->front(); }

      inline double w_max() const { return wvec_->back(); }

      inline double q_min() const { return qvec_->front(); }

      inline double q_max() const { return qvec_->back(); }

      inline const std::vector<double>& w_grid() const { return *wvec_; }

      inline const std::vector<double>& q_grid() const { return *qvec_; }

      inline const std::vector<NuclearResponses>& get_responses() const
        { return *responses_; }

      inline size_t response_index( size_t iw, size_t iq ) const {
        size_t num_q_points = qvec_->size();
        size_t r_idx = ( num_q_points * iw ) + iq;
        return r_idx;
      }

      NuclearResponses interpolate(double w, double q) const;
 
    protected:

      std::shared_ptr< std::vector<double> > wvec_;

      std::shared_ptr< std::vector<double> > qvec_;

      std::shared_ptr< std::vector<NuclearResponses> > responses_;

      bool get_neighbor_indices( const std::vector<double>& vec,
        double val, size_t& lower_index, size_t& upper_index ) const;
  };
 
}
 
inline std::ostream& operator<<(std::ostream& out,
  const marley::ResponseTable& rt)
{
  for ( size_t iw = 0u; iw < rt.w_grid().size(); ++iw ) {
    for ( size_t iq = 0u; iq < rt.q_grid().size(); ++iq ) {
      out << rt.w_grid().at( iw ) << "  " << rt.q_grid().at( iq ) << "  ";
      size_t ir = rt.response_index( iw, iq );
      out << rt.get_responses().at( ir ) << '\n';
    }
  }
  return out;
}