/**
 * @file
 * Some useful functions to perform simulations.
 */

#pragma once

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

namespace darkprop {

/**
 * Run a simulation and collect events that cross a sphere.
 */
template<typename Vector3, typename Value>
std::vector<Event<Vector3, Value>> simulate_cross_sphere(
        Particle<Vector3, Value>& p,
        Earth<Vector3, Value>& earth,
        Value Tcut,         /**< cutoff energy */
        const Vector3& r0,  /**< center of the sphere */
        Value R,            /**< radius */
        RandomNumber<Value>& rn)
{
    std::vector<Event<Vector3, Value>> events;
    unsigned long scatter_num = 0;
    unsigned long long_track = 1000000;
    auto prev_p = p.toEvent(false);
    do {
        prev_p = p.toEvent(false);
        earth.propagate(p, rn);
        sphere_cross(r0, R, prev_p, p, events);
        if (p.in_earth) {
            scatter(p, earth, rn);
            scatter_num++;
            if (scatter_num > long_track) {
                std::cerr << "scatter more than " << long_track << " times;"
                    << " depth = " << (constant::rEarth - p.r.norm())/unit::km
                    << " km" << std::endl;
                long_track += 1000000;
            }
            if (isnanf(p.r.norm())) {
                std::cout << "p.r.norm = "
                          << p.r.norm() / unit::km << std::endl;
                std::cout << "p.v.norm = "
                          << p.v.norm() * unit::sec / unit::km << std::endl;
                std::cout << "p.p3.norm = " << p.p3.norm() << std::endl;
                exit(1);
            }
            // if (p.T < Tcut || p.r.norm() < R) {
            if (p.T < Tcut) {
                p.in_earth = false;
            }
        }
    } while (p.in_earth);
    return events;
}


/**
 * Run a simulation and collect crossing events on several depths.
 */
template<typename Vector3, typename Value>
std::vector<std::vector<Event<Vector3, Value>>> simulate_cross_depth(
        Particle<Vector3, Value>& p,
        Earth<Vector3, Value>& earth,
        Value Tcut,                /**< cutoff energy */
        std::vector<Value> depth,  /**< collect events on these depth */
        RandomNumber<Value>& rn)
{
    std::vector<std::vector<Event<Vector3, Value>>> events(depth.size());
    unsigned long scatter_num = 0;
    unsigned long long_track = 1000000;
    Vector3 r0(0., 0., 0.);
    auto prev_p = p.toEvent(false);
    do {
        prev_p = p.toEvent(false);
        earth.propagate(p, rn);
        for (size_t i = 0; i != depth.size(); ++i) {
            sphere_cross(r0, constant::rEarth-depth[i], prev_p, p, events[i]);
        }
        if (p.in_earth) {
            scatter(p, earth, rn);
            scatter_num++;
            if (scatter_num > long_track) {
                std::cerr << std::scientific;
                std::cerr << "scatter more than " << long_track << " times;"
                    << " depth = " << (constant::rEarth - p.r.norm())/unit::km
                    << " km" << " T = " << p.T << " GeV" <<std::endl;
                long_track += 1000000;
            }
            if (p.T < Tcut) {
                p.in_earth = false;
            }
        }
    } while (p.in_earth);
    return events;
}


/**
 * Run a simulation and return the complete trajectory.
 */
template<typename Vector3, typename Value>
std::vector<Event<Vector3, Value>> simulate_track(
        Particle<Vector3, Value>& p,
        Earth<Vector3, Value>& earth,
        Value Tcut,  /**< cutoff energy */
        RandomNumber<Value>& rn)
{
    std::vector<Event<Vector3, Value>> events;
    unsigned long scatter_num = 0;
    unsigned long long_track = 1000000;
    do {
        events.push_back(p.toEvent());
        earth.propagate(p, rn);
        if (p.in_earth) {
            scatter(p, earth, rn);
            scatter_num++;
            if (scatter_num > long_track) {
                std::cerr << "scatter more than " << long_track << " times;"
                    << " depth = " << (constant::rEarth - p.r.norm())/unit::km
                    << " km" << std::endl;
                long_track += 1000000;
            }
            if (p.T < Tcut) {
                p.in_earth = false;
            }
        }
    } while (p.in_earth);
    events.push_back(p.toEvent());
    return events;
}

template<typename Vector3, typename Value>
std::vector<std::size_t> analyse_track(
        Particle<Vector3, Value>& p,
        Earth<Vector3, Value>& earth,
        Value Tcut,  /**< cutoff energy */
        Value depth,
        RandomNumber<Value>& rn)
{
    std::vector<Event<Vector3, Value>> events;
    std::size_t scatter_num = 0;
    std::size_t free_particles = 0;
    std::size_t below_vcut = 0;
    do {
        auto prev_p = p.toEvent(false);
        earth.propagate(p, rn);
        sphere_cross({0, 0, 0}, constant::rEarth - depth, prev_p, p, events);

        if (p.in_earth) {
            scatter(p, earth, rn);
            scatter_num++;
            if (p.T < Tcut) {
                p.in_earth = false;
                below_vcut++;
            }
        }
    } while (p.in_earth);
    std::size_t depth_crossing = events.size();
    if (scatter_num == 0) {
        free_particles++;
    }
    return std::vector<std::size_t> {below_vcut, free_particles,
                                     scatter_num, depth_crossing};
}

} // namespace darkprop
