#include <cmath>
#include <iostream>

#include <TApplication.h>
#include <TCanvas.h>
#include <TStyle.h>

#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Random.hh"
#include "Garfield/RandomEngineRoot.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewField.hh"

using namespace Garfield;

int main(int argc, char* argv[]) {
  TApplication app("app", &argc, argv);

  // Random seed for reproducibility.
 // RandomEngineRoot randomEngine(123456);
 // Random::SetEngine(randomEngine);

  // --------------------------------------------------
  // Gas
  // --------------------------------------------------
  MediumMagboltz gas("ar", 93., "co2", 5., "ic4h10", 2.);
  gas.SetTemperature(293.15);
  gas.SetPressure(760.);

  // --------------------------------------------------
  // Geometry: drift region + amplification region
  // --------------------------------------------------
  ComponentAnalyticField cmpD;
  ComponentAnalyticField cmpA;
  cmpD.SetMedium(&gas);
  cmpA.SetMedium(&gas);

  constexpr double bfield = 0.;
  cmpD.SetMagneticField(0., 0., bfield);
  cmpA.SetMagneticField(0., 0., bfield);

  // Amplification gap [cm]
  constexpr double yMesh = 0.0128;
  constexpr double vMesh = -500.;
  cmpA.AddPlaneY(yMesh, vMesh, "mesh");
  cmpA.AddPlaneY(0., 0., "bottom");

  // Drift gap [cm]
  constexpr double yDrift = yMesh + 0.7;
  constexpr double vDrift = -800.;
  cmpD.AddPlaneY(yDrift, vDrift, "top");
  cmpD.AddPlaneY(yMesh, vMesh, "mesh");



  // --------------------------------------------------
  // Sensor
  // --------------------------------------------------
  Sensor sensor;
  sensor.AddComponent(&cmpD);
  sensor.AddComponent(&cmpA);
  sensor.SetArea(-1, 0., -1., 1, yDrift, 1.);

  // --------------------------------------------------
  // Avalanche object for statistics
  // --------------------------------------------------
  AvalancheMicroscopic aval(&sensor);
  aval.EnableExcitationMarkers(false);
  aval.EnableIonisationMarkers(false);

  // --------------------------------------------------
  // Parameters of the primary electrons
  // --------------------------------------------------
  constexpr int nPrimaryElectrons = 5000;

  // Start slightly above the mesh, in the drift region.
  constexpr double dyStart = 1.e-4;   // [cm] = 1 micron above mesh
  const double yStart = yMesh + dyStart;

  // Fixed starting position.
  const double xStart = 0.;   // [cm]
  const double zStart = 0.0;    // [cm]
  const double tStart = 0.0;    // [ns]
  const double eStart = 1;    // [eV], small initial energy

  // Counters
  int nReachedMesh = 0;
  int nTransferred = 0;
  long long totalAvalancheElectrons = 0;
  long long totalAvalancheIons = 0;

  // --------------------------------------------------
  // Simulate 5000 primary electrons
  // --------------------------------------------------
  for (int i = 0; i < nPrimaryElectrons; ++i) {
    if (i % 100 == 0) {
      std::cout << i << " / " << nPrimaryElectrons << std::endl;
    }

    // Stage 1:
    // drift the electron from just above the mesh down to the mesh.
    aval.AvalancheElectron(xStart, yStart, zStart, tStart,
                           eStart, 0., -1., 0.);

    const auto& electrons1 = aval.GetElectrons();
    if (electrons1.empty()) continue;

    // Endpoint of the initial electron track.
    const auto& p1 = electrons1.front().path.back();

    // We only keep electrons that really reached the mesh plane.
    if (std::fabs(p1.y - yMesh) > 1.e-6) continue;
    ++nReachedMesh;

    // Stage 2:
    // move the electron slightly below the mesh and simulate avalanche.
    aval.AvalancheElectron(p1.x, yMesh - 1.e-6, p1.z, p1.t,
                           p1.energy, 0., -1., 0.);
    ++nTransferred;

    int ne = 0, ni = 0;
    aval.GetAvalancheSize(ne, ni);
    totalAvalancheElectrons += ne;
    totalAvalancheIons += ni;
  }

  // --------------------------------------------------
  // Print summary
  // --------------------------------------------------
  std::cout << "\n========== Summary ==========\n";
  std::cout << "Primary electrons started: " << nPrimaryElectrons << "\n";
  std::cout << "Reached mesh plane:        " << nReachedMesh << "\n";
  std::cout << "Transferred below mesh:    " << nTransferred << "\n";

  if (nTransferred > 0) {
    std::cout << "Mean avalanche electrons:  "
              << double(totalAvalancheElectrons) / nTransferred << "\n";
    std::cout << "Mean avalanche ions:       "
              << double(totalAvalancheIons) / nTransferred << "\n";
  }
  std::cout << "=============================\n";






  return 0;
}