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/Level.hh"
 
namespace marley {
 
  class MatrixElement {
 
    public:
 
      enum TransitionType {
        FERMI = 0, 
 
        GAMOW_TELLER = 1 
      };
 
      inline MatrixElement(double level_energy, double strength, TransitionType
        type, marley::Level* final_level = nullptr)
        : level_energy_(level_energy), strength_(strength), type_(type),
        final_level_(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 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_;
 
      TransitionType type_;
 
      marley::Level* final_level_;
  };
}