#include <cstdlib>
#include <fstream>
#include <iostream>
#include <math.h>
#include <sys/stat.h>

#include <TApplication.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <TROOT.h>

#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/ComponentNeBem3d.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/MediumConductor.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/SolidBox.hh"
#include "Garfield/SolidWire.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewGeometry.hh"
#include "Garfield/ViewMedium.hh"
#include "Garfield/ViewSignal.hh"

using namespace Garfield;

int main(int argc, char *argv[])
{

  TApplication app("app", &argc, argv);

  // load the gas medium for the internal space
  MediumMagboltz gas;
  gas.LoadGasFile("../gas_tables/ar_85_co2_15_1atm.gas");
  // load the Ar+ ion mobility files
  auto installdir = std::getenv("GARFIELD_INSTALL");
  if (!installdir) {
    std::cerr << "GARFIELD_INSTALL variable not set.\n";
    return 1;
  }
  const std::string path = installdir;
  gas.LoadIonMobility(path + "/share/Garfield/Data/IonMobility_Ar+_Ar.txt");

  std::cout << "Mass density: " << gas.GetMassDensity() << std::endl;

  // ViewMedium mediumView;
  // mediumView.SetMedium(&gas);
  // mediumView.SetRangeE(1e3,1e5,true);
  // mediumView.PlotElectronTownsend();

  // define the conductor material for the strips and wire
  MediumConductor plane_medium;
  MediumConductor strip_medium;
  MediumConductor wire_medium;

  // Defining the finite plane geometry through a GeometrySimple object
  GeometrySimple geo;

  // geometric parameters
  //----whole cell----
  // half-length of the assembly along z [cm]
  const double hlength = 100.;
  // gap half-length
  const double half_gap = 0.1;
  // half width of a cell
  const double cell_width = 1.;
  // half thickness of a cell
  const double cell_thick = 0.5 * cell_width;

  //----strips----
  // half dimensions
  const double stripWidth = 1 * cell_width - half_gap;
  const double stripThick = 35e-4;

  //----wires----
  const double diameter = 40;
  // signal wire radius [cm]
  const double r_signal_wire = (diameter * 1.e-4) / 2;
  // field wire radius
  const double r_field_wire = 100.e-4 / 2;

  //----Voltages----
  // strip voltage
  const double vStrip = 0.;
  // signal wire voltage
  const double v_signal_wire = +2000.;
  // field wire voltage
  const double v_field_wire = -2000.;

  // define the cathode strips
  SolidBox c_strip_up(0., cell_thick - stripThick, 0., stripWidth, stripThick,
                      hlength);
  SolidBox c_strip_down(0., -cell_thick + stripThick, 0., stripWidth,
                        stripThick, hlength);
  // define the signal wire
  SolidWire sens_wire(0., 0., 0., r_signal_wire, hlength);
  // set this wire as the electrode
  sens_wire.SetLabel("s");

  // define the field wires: left side
  SolidWire field_wire_u_l(-cell_width, 0., 0., r_field_wire, hlength);
  // define the field wires: right side
  SolidWire field_wire_u_r(+cell_width, 0., 0., r_field_wire, hlength);

  // set the boundary potentials for the components
  sens_wire.SetBoundaryPotential(v_signal_wire);
  field_wire_u_l.SetBoundaryPotential(v_field_wire);
  field_wire_u_r.SetBoundaryPotential(v_field_wire);
  c_strip_up.SetBoundaryPotential(vStrip);
  c_strip_down.SetBoundaryPotential(vStrip);

  // add wires strips and planes to the geometry
  geo.AddSolid(&sens_wire, &wire_medium);
  geo.AddSolid(&field_wire_u_l, &wire_medium);
  geo.AddSolid(&field_wire_u_r, &wire_medium);
  geo.AddSolid(&c_strip_up, &strip_medium);
  geo.AddSolid(&c_strip_down, &strip_medium);

  // add the boundary volume filled with gas
  geo.SetMedium(&gas);

  // set up the save folder
  const char *folder_name = Form(
      "../Debug/single_track_drift_simu/"
      "%2.0fumw_vSi%2.0f_vSt%2.0f_vF%2.0f",
      diameter, v_signal_wire, vStrip, v_field_wire);
  int check = mkdir(folder_name, 0777);
  if (!check)
    std::cout << "- - New folder created - -" << std::endl;
  else
    std::cout << "- - No new folder was created - -" << std::endl;

  // Before moving on one needs to compute the field with neBEM
  ComponentNeBem3d nebem;
  nebem.SetGeometry(&geo);
  nebem.SetTargetElementSize(0.1);
  nebem.SetMinMaxNumberOfElements(20, 20);
  nebem.SetNumberOfThreads(8);
  nebem.SetPeriodicityX(2 * cell_width);
  // nebem.SetPeriodicityY(2 * cell_thick);
  nebem.Initialise();

  // Use a Sensor object to interface to the transport classes
  Sensor sensor;
  // add the field and electrode components
  sensor.AddComponent(&nebem);
  sensor.AddElectrode(&nebem, "s");

  // setup the binning for the signal calculation (start, bin width[ns],bin
  // number)
  sensor.SetTimeWindow(0., 0.1, 3000);
  sensor.SetArea(-cell_width, -cell_thick, -hlength, cell_width, cell_thick,
                 hlength);

  // Simulating the ionization with Heed
  // first a track object is defined
  std::cout << "-- Particle track simulation --" << std::endl;
  TrackHeed track;
  track.SetParticle("muon");
  track.SetMomentum(2e9);
  track.SetSensor(&sensor);
  // compute the track
  track.EnableDebugging();

  // compute the electron drift lines
  DriftLineRKF drift;
  drift.SetSensor(&sensor);
  drift.UseWeightingPotential(true);
  drift.EnableAvalanche(true);
  drift.EnableIonTail(true);
  drift.EnableSignalCalculation();
  drift.EnableTownsendMap();

  // Setting up the viewer of the simulated drift lines
  // create a new canvas
  TCanvas *cD = new TCanvas("cD", " ", 1100, 500);
  // the ViewDrift object plots the cluster coordinates and drift line points
  ViewDrift driftView;
  // set the object with which to plot and the canvas
  driftView.SetCanvas(cD);
  // ava.EnablePlotting(&driftView);
  drift.EnablePlotting(&driftView);
  track.EnablePlotting(&driftView);

  // set the initial position: middle-left side of the cell
  const double x0 = -cell_width * 0.5, t0 = 0., y0 = -cell_thick, z0 = 0.;
  // set the track direction and define the track
  const double dx0 = 0., dy0 = 1., dz0 = 0.;

  track.NewTrack(x0, y0, z0, t0, dx0, dy0, dz0);
  // cluster is a Cluster struct instance: it has photons, ions and electrons as
  // attributes
  for (const auto &cluster : track.GetClusters()) {
    // declare endpoint coordinates
    double xe = 0., ye = 0., ze = 0., te = 0., ne = 0, ni = 0;
    double ex = 0., ey = 0., ez = 0.;
    Medium *m;
    int st = 0;
    // loop over the electrons in the cluster
    std::cout << "Cluster ne: " << (cluster.electrons).size()
              << ", nI: " << (cluster.ions).size() << std::endl;
    for (const auto &electron : cluster.electrons) {
      drift.DriftElectron(electron.x, electron.y, electron.z, electron.t);
      drift.GetEndPoint(xe, ye, ze, te, st);
      nebem.ElectricField(xe, ye, ze, ex, ey, ez, m, st);
      drift.GetAvalancheSize(ne, ni);
      std::cout << "end r: " << sqrt(xe * xe + ye * ye) << ", te: " << te
                << ", "
                << "e gain: " << drift.GetGain() << ", av_size: (" << ne << ", "
                << ni << "), E: " << sqrt(ex * ex + ey * ey + ez * ez)
                << std::endl;
    }
    for (const auto &ion : cluster.ions) {
      drift.DriftIon(ion.x, ion.y, ion.z, ion.t);
      drift.GetEndPoint(xe, ye, ze, te, st);
      std::cout << "Endpoint: (" << xe << ", " << ye << ", " << ze << ", " << te
                << "), "
                << "ion gain: " << drift.GetGain() << std::endl;
    }
  }
  //  plot the induced current on the wire by the electrons with ViewSignal
  ViewSignal signalView;
  signalView.SetSensor(&sensor);
  TCanvas *cD_signal = new TCanvas("cD_signal", " ", 900, 800);
  signalView.SetCanvas(cD_signal);
  signalView.PlotSignal("s", "ei");

  auto legend = new TLegend(0.7, 0.25, 0.8, 0.5);
  cD_signal->cd();
  legend->SetHeader("#splitline{Single track signal}{p=2 GeV #mu}", "C");
  legend->Draw();
  cD_signal->SaveAs(
      Form("%s/single_track_signal_e_ion_%2.0fumw.pdf", folder_name, diameter));

  std::cout << "Total induced charge: " << signalView.GetHistogram()->Integral()
            << std::endl;
  // set the object with which to plot and the canvas
  ViewGeometry geoView;
  geoView.SetGeometry(&geo);
  geoView.SetCanvas(cD);
  geoView.Plot2d();

  constexpr bool twod = true;
  constexpr bool drawaxis = false;
  driftView.Plot(twod, drawaxis);

  std::cout << "End of program" << std::endl;

  app.Run(true);
  return 0;
}