/* This file is part of MAUS: http://micewww.pp.rl.ac.uk/projects/maus
 *
 * MAUS is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * MAUS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with MAUS.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#ifndef recon_global_materialmodel_hh_
#define recon_global_materialmodel_hh_

#include <string>
#include <map>
#include <iostream>

class G4Material;

namespace MAUS {
class GeometryNavigator;

/** @class MaterialModel provides mean process model for a G4Material
 *
 *  MaterialModel calculates the mean energy loss, scatter and straggle for a
 *  given G4Material.
 *
 *  MaterialModel class holds sets of materials (stored by string name) that
 *  are explicitly enabled or disabled. Materials that are disabled will be
 *  ignored. The ignoring logic is handled by the children. By default all
 *  materials are disabled.
 */
class MaterialModel {
  public:
    /** Constructor (does nothing)
     */
    MaterialModel() {}
    /** Copy constructor not implemented */
    MaterialModel& operator=(const MaterialModel& mat);
    /** Destructor does nothing; note _material is a borrowed pointer */
    virtual ~MaterialModel() {}
    /** Inheritable copy constructor */
    virtual MaterialModel* Clone() = 0;

    /** Set the material to be modelled
     *
     *  @param material: the material to be modelled; it is a borrowed reference
     *                   (Geant4 usually will own this memory)
     */
    virtual void SetMaterial(const G4Material* material);

    /** Calculate the energy loss using Bethe Bloch formula
     *
     *  @param energy: energy of the incident particle
     *  @param mass: mass of the incident particle
     *  @param charge: charge of the incident particle
     *
     *  Takes the sum over elements in the target material; energy loss is given
     *  by the sum of energy loss due to each element according to Bethe Bloch.
     */
    virtual double dEdx(double energy, double mass, double charge);

    /** Calculate the energy growth due to curvature of Bethe Bloch
     *
     *  @param energy: energy of the incident particle
     *  @param mass: mass of the incident particle
     *  @param charge: charge of the incident particle
     *
     *  Performs numerical derivative of the Bethe Bloch curve by doing
     *  d2E/dEdz = [dE/dz(E+dE) - dE/dz(E-dE)]/2/dE with dE given by the
     *  d2EdxdE_delta_const.
     */
    virtual double d2EdxdE(double energy, double mass, double charge);

    /** Calculate the mean scattering angle squared per unit length.
     *
     *  @param energy: energy of the incident particle
     *  @param mass: mass of the incident particle
     *  @param charge: charge of the incident particle
     *
     *  Use theta_0 ~ const / beta / p sqrt(x/x0) so dtheta^2/dz ~ const / beta / p / x0
     */
    virtual double dtheta2dx(double energy, double mass, double charge);

    /** Calculate the mean energy straggling
     *
     *  @param energy: energy of the incident particle
     *  @param mass: mass of the incident particle
     *  @param charge: charge of the incident particle
     *
     *  Use estrag2 = const * gamma^2*(1-beta^2/2.)
     */
    virtual double estrag2(double energy, double mass, double charge);

    /** Set the material to be one corresponding to a given position x, y, z
     */
    virtual void SetMaterial(double x, double y, double z) = 0;

    /** Get the stored material
     *
     *  @returns the stored material; note this is a borrowed pointer. Geant4
     *           usually owns this memory.
     */
    virtual const G4Material* GetMaterial() const {return _material;}

    /** Return true if the material is enabled
     *
     *  @param material_name: the name of the material to be checked
     */
    static bool IsEnabledMaterial(std::string material_name);

    /** Enable a given material. Default is disabled.
     *
     *  @param material_name: the name of the material to be enabled
     */
    static void EnableMaterial(std::string material_name);

    /** Disable a given material.
     *
     *  @param material_name: the name of the material to be disabled
     */
    static void DisableMaterial(std::string material_name);

    /** Clear the enabled materials list */
    static void ClearEnabledMaterials();

    /** Set the default enablement (for items not in the list) */
    static void SetDefaultMaterialEnablement(bool default_enablement);

    /** Print a list of all enabled materials
     *
     *  @param out: prints the materials to the ostream out
     */
    static std::ostream& PrintMaterials(std::ostream& out);

    /** Set the energy loss calculation routine
     *  
     *  @param will_use_mean_energy_loss If true, use mean energy loss for dEdx 
     *  calculation; otherwise use most probable energy loss
     */
    void SetWillUseMeanEnergyLoss(bool will_use_mean_energy_loss);

    /** Get the energy loss calculation routine
     *  
     *  @returns will_use_mean_energy_loss If true, uses mean energy loss for
     *   dEdx calculation; otherwise uses most probable energy loss
     */
    bool GetWillUseMeanEnergyLoss() const;

    /** Set the material thickness; used for most probable energy loss routine
     *
     *  @param thickness material thickness in mm
     */
    void SetThickness(double thickness) {_thickness = thickness;}

    /** Get the material thickness; used for most probable energy loss routine
     *
     *  @returns material thickness in mm
     */
    double GetThickness() const {return _thickness;}
  protected:

    virtual double MeandEdx(double energy, double mass, double charge);
    virtual double MostProbabledEdx(double energy, double mass, double charge);
    // _material is a borrowed reference; data is owned by G4MaterialTable
    const G4Material* _material = NULL;
    GeometryNavigator* _navigator = NULL;
    double _n_e = 0.;
    double _I = 1.;
    double _x_0 = 0.;
    double _x_1 = 0.;
    double _C = 0.;
    double _a = 0.;
    double _k = 0.;
    double _rad_len_ratio = 0.;
    double _density = 0.;
    double _z_over_a = 0.;
    double _thickness = 1.;

    // if true, use mean energy loss for dEdx calculation; otherwise use
    // most probable energy loss
    bool _will_use_mean_energy_loss = true;
    // map of enabled materials; if not in the map, default enablement is used
    static std::map<std::string, bool> _enabled_materials;
    static bool _default_enablement;

    static constexpr double _d2EdxdE_delta_const = 0.1;
    // 2.1047291091867513e-13; // 4 pi (r_e m_e c^2)^2
    static constexpr double _estrag_const = 0.157;
    static constexpr double _dedx_constant = 0.307075;
    static constexpr double _m_e = 0.510998928;
    static constexpr double c_l = 299.792458; // mm*ns^{-1}
};
}

#endif