#include <iostream>
#include <iomanip>

#include "Garfield/Sensor.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/TrackSrim.hh"
#include "Garfield/ComponentConstant.hh"

using namespace Garfield;

int main() {
  // ------------------------------------------------------------
  // Gas: He/N2 mixture at 6 atm.
  // ------------------------------------------------------------
  MediumMagboltz gas;
  const std::string gasfile = "../gas/he_98_n2_2_6atm_penning.gas";

  if (!gas.LoadGasFile(gasfile)) {
    std::cerr << "Failed to load gas file: " << gasfile << "\n";
    return 1;
  }

  // ------------------------------------------------------------
  // Minimal constant-field component.
  // This provides a valid medium for TrackSrim.
  // The electric field is not important here because we only test
  // SRIM cluster generation, not charge drift.
  // ------------------------------------------------------------
  ComponentConstant cmp;
  cmp.SetArea(-1.0, -1.0, 0.0, 1.0, 1.0, 2.0); // cm
  cmp.SetMedium(&gas);
  cmp.SetElectricField(0.0, 0.0, 0.0);          // V/cm

  Sensor sensor;
  sensor.AddComponent(&cmp);
  sensor.SetArea(-1.0, -1.0, 0.0, 1.0, 1.0, 2.0); // cm

  // ------------------------------------------------------------
  // Load SRIM table for Li in gas.
  // ------------------------------------------------------------
  TrackSrim trLi(&sensor);

  const std::string srimFile =
      "../trim/Lithium_in_Helium_6T_range (gas).txt";

  if (!trLi.ReadFile(srimFile)) {
    std::cerr << "Failed to read SRIM file: " << srimFile << "\n";
    return 1;
  }

  trLi.SetTargetClusterSize(1);
  trLi.SetClustersMaximum(10000);

  // ------------------------------------------------------------
  // Test Li7 track.
  // Energy is given in eV.
  // Position is in cm.
  // ------------------------------------------------------------
  const double x0 = 0.0;
  const double y0 = 0.0;
  const double z0 = 0.6001; // cm, inside [0, 2] cm

  const double t0 = 0.0; // ns

  const double dx = 0.0;
  const double dy = 0.0;
  const double dz = 1.0;

  const double kineticEnergyEv = 20.0e3; // 20 keV Li7

  trLi.SetKineticEnergy(kineticEnergyEv);

  if (!trLi.NewTrack(x0, y0, z0, t0, dx, dy, dz)) {
    std::cerr << "TrackSrim::NewTrack failed.\n";
    return 1;
  }

  const auto& clusters = trLi.GetClusters();

  if (clusters.empty()) {
    std::cout << "No clusters were generated.\n";
    return 0;
  }

  double totalClusterEnergyEv = 0.0;
  unsigned long long totalElectrons = 0;

  std::cout << std::scientific << std::setprecision(6);

  std::cout << "Number of clusters: " << clusters.size() << "\n";

  for (size_t i = 0; i < clusters.size(); ++i) {
    const auto& c = clusters[i];

    totalClusterEnergyEv += c.energy;
    totalElectrons += c.n;

    std::cout << "Cluster " << i
              << ": z = " << c.z << " cm"
              << ", cluster.energy = " << c.energy << " eV"
              << ", cluster.kinetic = " << c.kinetic << " MeV"
              << ", n = " << c.n
              << "\n";
  }

  const auto& last = clusters.back();

  std::cout << "\nSummary\n";

  std::cout << "Initial Li7 kinetic energy: "
            << kineticEnergyEv << " eV\n";

  std::cout << "Total energy converted into clusters: "
            << totalClusterEnergyEv << " eV\n";

  std::cout << "Total number of electron-ion pairs: "
            << totalElectrons << "\n";

  if (totalElectrons > 0) {
    std::cout << "Mean cluster energy per electron-ion pair: "
              << totalClusterEnergyEv / totalElectrons
              << " eV/pair\n";
  }

  std::cout << "Last cluster position z: "
            << last.z << " cm\n";

  std::cout << "Last cluster kinetic energy: "
            << last.kinetic * 1.0e6 << " eV\n";

  std::cout << "Last cluster kinetic energy minus last cluster.energy: "
            << last.kinetic * 1.0e6 - last.energy
            << " eV\n";

  return 0;
}