MatrixElement.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
 
#include "marley/Error.hh"
#include "marley/Level.hh"
 
namespace marley {

  class MatrixElement {
 
    public:

      enum TransitionType {
        FERMI = 0, 
 
        GAMOW_TELLER = 1 
      };

      inline MatrixElement(double level_energy, double strength, TransitionType
        type, double err_low, double err_high,
        marley::Level* final_level = nullptr)
        : level_energy_(level_energy), strength_(strength),
        strength_err_low_(err_low), strength_err_high_(err_high),
        type_(type), final_level_(final_level) {}

      inline MatrixElement(double level_energy, double strength, TransitionType
        type, marley::Level* final_level = nullptr)
        : MatrixElement(level_energy, strength, type, 0., 0., final_level) {}

      inline double level_energy() const {
        // If this matrix element accesses a discrete final nuclear level,
        // then return that level's excitation energy
        if ( final_level_ ) return final_level_->energy();
        // Otherwise, we're in the ubound continuum, and we can just
        // use the tabulated value from the reaction matrix element data file
        else return level_energy_;
      }

      inline double tabulated_level_energy() const { return level_energy_; }

      inline double strength() const { return strength_; }

      inline double strength_err_low() const { return strength_err_low_; }

      inline double strength_err_high() const { return strength_err_high_; }

      inline TransitionType type() const { return type_; }

      inline const marley::Level* level() const { return final_level_; }

      inline marley::Level* level() { return final_level_; }

      inline void set_level_energy(double energy) { level_energy_ = energy; }

      inline void set_strength(double strength) { strength_ = strength; }

      inline void set_type(TransitionType type) { type_ = type; }

      inline void set_level(marley::Level* lev) { final_level_ = lev; }

      inline double cos_theta_pdf(double cos_theta_c_cm, double beta_c_cm) const
      {
        double pdf = 0.;
        if ( type_ == TransitionType::FERMI ) {
          pdf = 1. + beta_c_cm * cos_theta_c_cm;
        }
        else if ( type_ == TransitionType::GAMOW_TELLER ) {
          pdf = (3. - beta_c_cm * cos_theta_c_cm) / 3.;
        }
        else throw marley::Error("Unrecognized transition type encountered"
          " in marley::MatrixElement::cos_theta_pdf()");
 
        // Normalize to unit integral
        pdf *= 0.5;
        return pdf;
      }

      inline std::string type_str() const {
        if ( type_ == TransitionType::FERMI ) return "Fermi";
        else if ( type_ == TransitionType::GAMOW_TELLER ) return "Gamow-Teller";
        else throw marley::Error( "Unrecognized transition type encountered"
          " in marley::MatrixElement::type_str()" );
      }
 
    protected:

      double level_energy_;

      double strength_;

      double strength_err_low_ = 0.;

      double strength_err_high_ = 0.;

      TransitionType type_;

      marley::Level* final_level_;
  };
}