src/physics/ImpactIonization.hpp

#ifndef ALLPIX_IMPACTIONIZATION_MODELS_H
#define ALLPIX_IMPACTIONIZATION_MODELS_H

#include <limits>
#include <typeindex>

#include <TFormula.h>

#include "exceptions.h"

#include "core/config/Configuration.hpp"
#include "core/utils/log.h"
#include "core/utils/unit.h"
#include "objects/SensorCharge.hpp"

namespace allpix {

    class ImpactIonizationModel {
    public:
        explicit ImpactIonizationModel(double threshold) : threshold_(threshold){};

        virtual ~ImpactIonizationModel() = default;

        virtual double operator()(const CarrierType& type, double efield_mag, double step) const {
            if(std::fabs(efield_mag) < threshold_) {
                return 1.;
            }
            return std::exp(step * gain_factor(type, efield_mag));
        };

    protected:
        virtual double gain_factor(const CarrierType& type, double efield_mag) const = 0;
        double threshold_{std::numeric_limits<double>::max()};
    };

    class NoImpactIonization : virtual public ImpactIonizationModel {
    public:
        NoImpactIonization() : ImpactIonizationModel(std::numeric_limits<double>::max()){};
        double operator()(const CarrierType&, double, double) const override { return 1.; };

    private:
        double gain_factor(const CarrierType&, double) const override { return 1.; };
    };

    class Massey : virtual public ImpactIonizationModel {
    public:
        Massey(double temperature, double threshold)
            : ImpactIonizationModel(threshold), electron_a_(Units::get(4.43e5, "/cm")),
              electron_b_(Units::get(9.66e5, "V/cm") + Units::get(4.99e2, "V/cm/K") * temperature),
              hole_a_(Units::get(1.13e6, "/cm")),
              hole_b_(Units::get(1.71e6, "V/cm") + Units::get(1.09e3, "V/cm/K") * temperature) {}

    private:
        double gain_factor(const CarrierType& type, double efield_mag) const override {
            if(type == CarrierType::ELECTRON) {
                return electron_a_ * std::exp(-1. * electron_b_ / efield_mag);
            } else {
                return hole_a_ * std::exp(-1. * hole_b_ / efield_mag);
            }
        };

    protected:
        double electron_a_;
        double electron_b_;

        double hole_a_;
        double hole_b_;
    };

    class MasseyOptimized : virtual public Massey {
    public:
        MasseyOptimized(double temperature, double threshold)
            : ImpactIonizationModel(threshold), Massey(temperature, threshold) {
            electron_a_ = Units::get(1.186e6, "/cm");
            electron_b_ = Units::get(1.020e6, "V/cm") + Units::get(1.043e3, "V/cm/K") * temperature;
            hole_a_ = Units::get(2.250e6, "/cm");
            hole_b_ = Units::get(1.851e6, "V/cm") + Units::get(1.828e3, "V/cm/K") * temperature;
        };
    };

    class VanOverstraetenDeMan : virtual public ImpactIonizationModel {
    public:
        VanOverstraetenDeMan(double temperature, double threshold)
            : ImpactIonizationModel(threshold),
              gamma_(std::tanh(Units::get(0.063, "eV") / (2. * Units::get(8.6173333e-5, "eV/K") * 300.)) /
                     std::tanh(Units::get(0.063, "eV") / (2. * Units::get(8.6173333e-5, "eV/K") * temperature))),
              e_zero_(Units::get(4.0e5, "V/cm")), electron_a_(Units::get(7.03e5, "/cm")),
              electron_b_(Units::get(1.231e6, "V/cm")), hole_a_low_(Units::get(1.582e6, "/cm")),
              hole_a_high_(Units::get(6.71e5, "/cm")), hole_b_low_(Units::get(2.036e6, "V/cm")),
              hole_b_high_(Units::get(1.693e6, "V/cm")) {}

    private:
        double gain_factor(const CarrierType& type, double efield_mag) const override {
            if(type == CarrierType::ELECTRON) {
                return gamma_ * electron_a_ * std::exp(-(gamma_ * electron_b_ / efield_mag));
            } else {
                auto a = (std::abs(efield_mag) > e_zero_ ? hole_a_high_ : hole_a_low_);
                auto b = (std::abs(efield_mag) > e_zero_ ? hole_b_high_ : hole_b_low_);
                return gamma_ * a * std::exp(-(gamma_ * b / efield_mag));
            }
        };

    protected:
        double gamma_;
        double e_zero_;

        double electron_a_;
        double electron_b_;

        double hole_a_low_;
        double hole_a_high_;
        double hole_b_low_;
        double hole_b_high_;
    };

    class VanOverstraetenDeManOptimized : virtual public VanOverstraetenDeMan {
    public:
        VanOverstraetenDeManOptimized(double temperature, double threshold)
            : ImpactIonizationModel(threshold), VanOverstraetenDeMan(temperature, threshold) {
            gamma_ = std::tanh(Units::get(0.0758, "eV") / (2. * Units::get(8.6173333e-5, "eV/K") * 300.)) /
                     std::tanh(Units::get(0.0758, "eV") / (2. * Units::get(8.6173333e-5, "eV/K") * temperature));
            electron_a_ = Units::get(1.149e6, "/cm");
            electron_b_ = Units::get(1.325e6, "V/cm");
            hole_a_low_ = Units::get(2.519e6, "/cm");
            hole_b_low_ = Units::get(2.428e6, "V/cm");
            // The publication uses the same parameters for low and high electric field regions:
            hole_a_high_ = Units::get(2.519e6, "/cm");
            hole_b_high_ = Units::get(2.428e6, "V/cm");
        };
    };

    class OkutoCrowell : virtual public ImpactIonizationModel {
    public:
        OkutoCrowell(double temperature, double threshold)
            : ImpactIonizationModel(threshold), electron_ac_(Units::get(0.426, "/V") * (1. + 3.05e-4 * (temperature - 300))),
              electron_bd_(Units::get(4.81e5, "V/cm") * (1. + 6.86e-4 * (temperature - 300))),
              hole_ac_(Units::get(0.243, "/V") * (1. + 5.35e-4 * (temperature - 300))),
              hole_bd_(Units::get(6.53e5, "V/cm") * (1. + 5.67e-4 * (temperature - 300))) {}

    private:
        double gain_factor(const CarrierType& type, double efield_mag) const override {
            if(type == CarrierType::ELECTRON) {
                return electron_ac_ * efield_mag * std::exp(-1 * electron_bd_ * electron_bd_ / efield_mag / efield_mag);
            } else {
                return hole_ac_ * efield_mag * std::exp(-1 * hole_bd_ * hole_bd_ / efield_mag / efield_mag);
            }
        };

    protected:
        double electron_ac_;
        double electron_bd_;
        double hole_ac_;
        double hole_bd_;
    };

    class OkutoCrowellOptimized : virtual public OkutoCrowell {
    public:
        OkutoCrowellOptimized(double temperature, double threshold)
            : ImpactIonizationModel(threshold), OkutoCrowell(temperature, threshold) {
            electron_ac_ = Units::get(0.289, "/V") * (1. + 9.03e-4 * (temperature - 300));
            electron_bd_ = Units::get(4.01e5, "V/cm") * (1. + 1.11e-3 * (temperature - 300));
            hole_ac_ = Units::get(0.202, "/V") * (1. - 2.20e-3 * (temperature - 300));
            hole_bd_ = Units::get(6.40e5, "V/cm") * (1. + 8.25e-4 * (temperature - 300));
        };
    };

    class Bologna : virtual public ImpactIonizationModel {
    public:
        Bologna(double temperature, double threshold)
            : ImpactIonizationModel(threshold),
              electron_a_(Units::get(4.3383, "V") + Units::get(-2.42e-12, "V") * std::pow(temperature, 4.1233)),
              electron_b_(Units::get(0.235, "V")),
              electron_c_(Units::get(1.6831e4, "V/cm") + Units::get(4.3796, "V/cm") * temperature +
                          Units::get(0.13005, "V/cm") * std::pow(temperature, 2)),
              electron_d_(Units::get(1.2337e6, "V/cm") + Units::get(1.2039e3, "V/cm") * temperature +
                          Units::get(0.56703, "V/cm") * std::pow(temperature, 2)),
              hole_a_(Units::get(2.376, "V") + Units::get(1.033e-2, "V") * temperature),
              hole_b_(Units::get(0.17714, "V") * std::exp(Units::get(-2.178e-3, "/K") * temperature)),
              hole_c_(Units::get(9.47e-3, "V/cm") * std::pow(temperature, 2.4924)),
              hole_d_(Units::get(1.4043e6, "V/cm") + Units::get(2.9744e3, "V/cm") * temperature +
                      Units::get(1.4829, "V/cm") * std::pow(temperature, 2)) {}

    private:
        double gain_factor(const CarrierType& type, double efield_mag) const override {
            if(type == CarrierType::ELECTRON) {
                return efield_mag / (electron_a_ + electron_b_ * std::exp(electron_d_ / (efield_mag + electron_c_)));
            } else {
                return efield_mag / (hole_a_ + hole_b_ * std::exp(hole_d_ / (efield_mag + hole_c_)));
            }
        };

        double electron_a_;
        double electron_b_;
        double electron_c_;
        double electron_d_;
        double hole_a_;
        double hole_b_;
        double hole_c_;
        double hole_d_;
    };

    class CustomGain : public ImpactIonizationModel {
    public:
        CustomGain(const Configuration& config, double threshold) : ImpactIonizationModel(threshold) {
            electron_gain_ = configure_gain(config, CarrierType::ELECTRON);
            hole_gain_ = configure_gain(config, CarrierType::HOLE);
        };

        double gain_factor(const CarrierType& type, double efield_mag) const override {
            if(type == CarrierType::ELECTRON) {
                return electron_gain_->Eval(efield_mag);
            } else {
                return hole_gain_->Eval(efield_mag);
            }
        };

    private:
        std::unique_ptr<TFormula> electron_gain_;
        std::unique_ptr<TFormula> hole_gain_;

        std::unique_ptr<TFormula> configure_gain(const Configuration& config, const CarrierType type) {
            std::string name = (type == CarrierType::ELECTRON ? "electrons" : "holes");
            auto function = config.get<std::string>("multiplication_function_" + name);
            auto parameters = config.getArray<double>("multiplication_parameters_" + name, {});

            auto gain = std::make_unique<TFormula>(("multiplication_" + name).c_str(), function.c_str());

            if(!gain->IsValid()) {
                throw InvalidValueError(config,
                                        "multiplication_function_" + name,
                                        "The provided model is not a valid ROOT::TFormula expression");
            }

            // Check if number of parameters match up
            if(static_cast<size_t>(gain->GetNpar()) != parameters.size()) {
                throw InvalidValueError(config,
                                        "multiplication_parameters_" + name,
                                        "The number of provided parameters and parameters in the function do not match");
            }

            // Set the parameters
            for(size_t n = 0; n < parameters.size(); ++n) {
                gain->SetParameter(static_cast<int>(n), parameters[n]);
            }

            return gain;
        };
    };

    class ImpactIonization {
    public:
        ImpactIonization() = default;

        explicit ImpactIonization(const Configuration& config) {
            try {
                auto model = config.get<std::string>("multiplication_model", "none");
                std::transform(model.begin(), model.end(), model.begin(), ::tolower);
                auto temperature = config.get<double>("temperature");
                auto threshold = config.get<double>("multiplication_threshold");

                if(model == "massey") {
                    model_ = std::make_unique<Massey>(temperature, threshold);
                } else if(model == "massey_optimized") {
                    model_ = std::make_unique<MasseyOptimized>(temperature, threshold);
                } else if(model == "overstraeten") {
                    model_ = std::make_unique<VanOverstraetenDeMan>(temperature, threshold);
                } else if(model == "overstraeten_optimized") {
                    model_ = std::make_unique<VanOverstraetenDeManOptimized>(temperature, threshold);
                } else if(model == "okuto") {
                    model_ = std::make_unique<OkutoCrowell>(temperature, threshold);
                } else if(model == "okuto_optimized") {
                    model_ = std::make_unique<OkutoCrowellOptimized>(temperature, threshold);
                } else if(model == "bologna") {
                    model_ = std::make_unique<Bologna>(temperature, threshold);
                } else if(model == "none") {
                    LOG(INFO) << "No impact ionization model chosen, charge multiplication not simulated";
                    model_ = std::make_unique<NoImpactIonization>();
                } else if(model == "custom") {
                    model_ = std::make_unique<CustomGain>(config, threshold);
                } else {
                    throw InvalidModelError(model);
                }
                LOG(INFO) << "Selected impact ionization model \"" << model << "\"";
            } catch(const ModelError& e) {
                throw InvalidValueError(config, "multiplication_model", e.what());
            }
        }

        template <class... ARGS> double operator()(ARGS&&... args) const {
            return model_->operator()(std::forward<ARGS>(args)...);
        }

        template <class T> bool is() const { return dynamic_cast<T*>(model_.get()) != nullptr; }

    private:
        std::unique_ptr<ImpactIonizationModel> model_{};
    };

} // namespace allpix

#endif