/**
 * @file
 * Spin independent elastic scattering DM model with constant cross section.
 */

#pragma once

#include <iostream>
#include <ios>
#include "Base.hpp"
#include "Const.hpp"
#include "Vector3d.hpp"

namespace darkprop {

/**
 * Spin independent elastic scattering DM model with constant cross section.
 */
template<typename Vector3=Vector3d<double>, typename Value=double>
class SIDM : public Particle<Vector3, Value>
{
public:
    // TODO Base class constructor
    SIDM(Value t_sigma0=1e-30*unit::cm*unit::cm,
         Value t_m=1.*unit::GeV,
         Value t_t=0.0) {
        this->t = t_t;
        this->sigma0 = t_sigma0;
        this->m = t_m;
        this->in_earth = false;
    }
    SIDM(Value t_sigma0,
         Value t_m,
         Value t_t,
         Vector3 t_r,
         Vector3 t_v)
    {
        this->sigma0 = t_sigma0;
        this->m = t_m;
        this->t = t_t;
        this->r = t_r;
        this->v = t_v;
        this->updateP3TwithV();
        this->updateEp();
        this->updateEr();
        this->in_earth = false;
    }
    SIDM(Vector3 t_r,
         Vector3 t_p,
         Value t_sigma0,
         Value t_m,
         Value t_t)
    {
        this->sigma0 = t_sigma0;
        this->m = t_m;
        this->t = t_t;
        this->r = t_r;
        this->p3 = t_p;
        this->T = this->TFromP3();
        this->updateVwithP3T();
        this->updateEp();
        this->updateEr();
        this->in_earth = false;
    }
    ~SIDM() {}
    Value sigmaEff(const Target& t) const {
        if (t.A == -1) {
            // for electron
            return this->sigma0;
        } else {
            return this->sigma0 * t.A*t.A * std::pow(t.m / constant::mu *
                        (this->m + constant::mu) / (this->m + t.m), 2.0);
        }
    }
    Value maximumRecoilT(const Target& target) const {
        return this->T * (this->T + 2. * this->m)
               / (this->T + std::pow(this->m + target.m, 2) / (2. * target.m));
    }
    Value differentialCrossSection(Value Tr, const Target& t) const {
        Value Trmax = maximumRecoilT(t);
        if (Tr > Trmax || Tr < 0) return 0.0;
        return sigmaEff(t) / Trmax;
    }
    Value totalCrossSection(const Target& t) const override {
        return sigmaEff(t);
    }
    // Value sampleTanTheta(const Target& t, RandomNumber<Value>& rn) {
    //     Value xi = rn.uniform_xi();

    //     Value coschi = 1. - 2.*xi;
    //     Value sinchi = std::sqrt(1. - coschi*coschi);
    //     Value m1plusm2 = this->m + t.m;
    //     Value E1plusm2 = this->T + m1plusm2;

    //     return sinchi
    //         / (E1plusm2/std::sqrt(2.*t.m*this->T+m1plusm2*m1plusm2)) //gamma
    //         / (coschi + (this->m*m1plusm2+this->T*t.m)/(t.m * E1plusm2));
    // }

    Value scatteringSampleTr(const Target& t, RandomNumber<Value>& rn) override
    {
        return maximumRecoilT(t) * rn.uniform_xi();
    }

    Value scatteringCosTheta(const Target& t, Value Tr) override {
        // it seems no way to ensure |costheta| <= 1 :(
        Value p12 = this->T * (this->T + 2. * this->m);  // p1 squared
        Value Tf = this->T - Tr;  // final energy of DM
        // Value costheta = (p12 - Tr * (this->T + this->m + t.m))
        //                  / std::sqrt(p12 * Tf * (Tf + 2. * this->m));

        Value numerator = p12 - Tr * (this->T + this->m + t.m);
        Value costheta = std::sqrt(std::pow(numerator, 2.0)
                         / p12 / (Tf * (Tf + 2. * this->m)));
        costheta = costheta * (numerator > 0 ? 1 : -1);

        if (costheta > 1.0) {
            if (costheta > 1.0 + 1e-15) {
                std::cerr << std::scientific;
                std::cerr.precision(18);
                std::cerr << "!!!Warning!!! costheta = " << costheta
                    << " mchi = " << this->m
                    << " T = " << this->T
                    << " Tr = " << Tr
                    << " Tf = " << Tf << std::endl;
                std::cerr << "set costheta = 1\n";
            }
            costheta = 1.0;
        } else if (costheta < -1.0) {
            if (costheta < -(1.0 + 1e-15)) {
                std::cerr << std::scientific;
                std::cerr.precision(18);
                std::cerr << "!!!Warning!!! costheta = " << costheta
                    << " mchi = " << this->m
                    << " T = " << this->T
                    << " Tr = " << Tr
                    << " Tf = " << Tf << std::endl;
                std::cerr << "set costheta = -1\n";
            }
            costheta = -1.0;
        }
        return costheta;
    }
};

} // namespace darkprop
