Reaction.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 <memory>
#include <string>
#include <vector>
 
#include "marley/CoulombCorrector.hh"
#include "marley/MassTable.hh"
#include "marley/TargetAtom.hh"
 
namespace HepMC3 {
  class GenEvent;
  class GenParticle;
}
 
namespace marley {
 
  class Generator;
  class JSON;
  class Parity;
  class StructureDatabase;

  class Reaction {
 
    public:
 
      virtual ~Reaction() = default;

      Reaction( const std::string& source_file );

      enum ProcessType {
        Unknown = -1, 
        NeutrinoCC_Discrete = 0,     
        AntiNeutrinoCC_Discrete = 1, 
        NC_Discrete = 2,             
        NuElectronElastic = 3, 
        NeutrinoCC_Continuum = 4,     
        AntiNeutrinoCC_Continuum = 5, 
        NC_Continuum = 6,             
        StandaloneDecay = 7, 
      };

      enum DataFormat {
        DiscreteStrengths = 0,
        MultipoleResponses = 1
      };

      virtual double total_xs( int pdg_a, double KEa ) const = 0;

      virtual std::shared_ptr< HepMC3::GenEvent > create_event( int pdg_a,
        double KEa, marley::Generator& gen ) const = 0;

      inline const std::string& get_description() const { return description_; }

      inline ProcessType process_type() const { return process_type_; }

      inline const std::string& source_file() const { return source_file_; }
 
      static std::string proc_type_to_string( const ProcessType& pt );

      inline int pdg_a() const { return pdg_a_; }

      inline int pdg_b() const { return pdg_b_; }

      virtual double threshold_kinetic_energy() const = 0;

      virtual marley::TargetAtom atomic_target() const = 0;

      static std::vector< std::unique_ptr< Reaction > >
        load_from_file( const std::string& filename, StructureDatabase& db,
        CoulombCorrector::CoulombMode coulomb_mode, const JSON& ff_config );

      static int get_ejectile_pdg( int pdg_a, ProcessType proc_type );

      static void get_residue_pdg_and_charge( ProcessType proc_type,
        int pdg_b, int& pdg_d, int& q_d );
 
    protected:
 
      int pdg_a_; 
      int pdg_b_; 
      int pdg_c_; 
      int pdg_d_; 
 
      double ma_; 
      double mb_; 
      double mc_; 
 
      // Mutable because, for the nuclear reaction case, we
      // may modify the final nuclear mass based on a sampled
      // excitation energy
      mutable double md_; 

      std::string description_;

      ProcessType process_type_;

      void two_two_scatter( double KEa, double& s, double& Ec_cm,
        double& pc_cm, double& Ed_cm ) const;

      virtual std::shared_ptr< HepMC3::GenEvent > make_event_object(
        double KEa, double pc_cm, double cos_theta_c_cm, double phi_c_cm,
        double Ec_cm, double Ed_cm, int residue_status ) const;

      virtual std::shared_ptr< HepMC3::GenEvent > make_event_object(
        double KEa, const std::shared_ptr< HepMC3::GenParticle >& ejectile,
        const std::shared_ptr< HepMC3::GenParticle >& residue ) const;

      static const std::vector<int>& get_projectiles( ProcessType proc_type );

      const std::string source_file_;
  };
 
}