#include <cstdlib>
#include <iostream>
#include <TROOT.h>
#include <TApplication.h>

#include "Garfield/MediumMagboltz.hh"
#include "Garfield/MediumConductor.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/TrackSrim.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/ViewMedium.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ViewCell.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/ViewGeometry.hh"
#include "Garfield/FundamentalConstants.hh"
#include "Garfield/Random.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/SolidBox.hh"
#include "Garfield/SolidWire.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/ComponentNeBem3d.hh"

using namespace Garfield;
using namespace std;

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

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

  // Setup the gas.
  MediumMagboltz gas;
  // Set the temperature [K] and pressure [Torr].
  gas.SetTemperature(288.15);
  gas.SetPressure(45.);
  // Set the gas mixture.
  // gas.SetFieldGrid(100., 100.e3, 20, true);
  gas.SetComposition("isobutane", 100);
  // const int ncoll = 10;
  // gas.GenerateGasTable(ncoll);
  // gas.WriteGasFile("isobutane.gas");
  gas.LoadGasFile("isobutane.gas");
  const std::string path = std::getenv("GARFIELD_INSTALL");
  gas.LoadIonMobility(path + "/share/Garfield/Data/IonMobility_C8Hn+_iC4H10.txt");
  gas.SetMaxElectronEnergy(200.);
  gas.EnableDrift();
  gas.EnablePrimaryIonisation();
  gas.Initialise();

  MediumConductor metal;

  // ViewMedium mediumView;
  // mediumView.SetMedium(&gas);
  // mediumView.PlotElectronVelocity('e');

  GeometrySimple geo;
  // Lenghts in [cm]
  constexpr double rWire = 10.e-4;
  constexpr double halfLA = 0.5 / 2;
  constexpr double halfLX = 0.8 / 2;
  constexpr double halfLY = 0.5 / 2;
  constexpr double dy = 0.24;
  // Wire spacings [cm]
  constexpr double periodA = 0.05;
  constexpr double periodX = 0.1;
  constexpr double periodY = 0.1;
  // Number of wires
  const int NWireX = 6;  //6
  const int NWireY = 9;  //9
  const int NWireA = 17; //17
  SolidWire *WireA[NWireA];
  SolidWire *WireX[NWireX];
  SolidWire *WireY[NWireY];
  // Anode voltage [V]
  constexpr double vAnode = 460.; //460.;

  for (int i = 0; i < NWireA; i++)
  {
    WireA[i] = new SolidWire(i * periodA, 0, halfLA, rWire, halfLA, 0, 0, 1);
    WireA[i]->SetBoundaryPotential(vAnode);
    // WireA[i]->SetBoundaryDielectric();
    geo.AddSolid(WireA[i], &metal);
  }
  for (int i = 0; i < NWireX; i++)
  {
    WireX[i] = new SolidWire(halfLX, dy, i * periodX, rWire, halfLX, 1, 0, 0);
    WireX[i]->SetLabel(Form("X%d", i));
    WireX[i]->SetBoundaryPotential(0.);
    // WireX[i]->SetBoundaryDielectric();
    geo.AddSolid(WireX[i], &metal);
  }
  for (int i = 0; i < NWireY; i++)
  {
    WireY[i] = new SolidWire(i * periodY, -1. * dy, halfLY, rWire, halfLY, 0, 0, 1);
    WireY[i]->SetLabel(Form("Y%d", i));
    WireY[i]->SetBoundaryPotential(0.);
    // WireY[i]->SetBoundaryDielectric();
    geo.AddSolid(WireY[i], &metal);
  }

  geo.SetMedium(&gas);

  ViewGeometry geoView;
  geoView.SetGeometry(&geo);
  geoView.Plot3d();

  // Construct electric field
  ComponentNeBem3d nebem;
  nebem.SetStoreReadOptions(1, 0, 1, 0, 1, 0, 1, 0, 1, 0);
  // nebem.SetStoreReadOptions(0, 1, 0, 1, 0, 1, 0, 1, 1, 0);
  nebem.SetGeometry(&geo);
  nebem.SetTargetElementSize(0.1);
  nebem.SetNumberOfThreads(10);
  // nebem.EnableDebugging();
  nebem.UseSVDInversion();
  // nebem.SetReuseModel();
  nebem.Initialise();

  // Medium *medium = nullptr;
  // double ex = 0., ey = 0., ez = 0., v = 0.;
  // int status = 0;
  // nebem.ElectricField(0, 0, 0, ex, ey, ez, v, medium, status);
  // std::printf("E = (%15.8f, %15.8f %15.8f), V = %15.8f, status = %d\n", ex, ey, ez, v, status);

  // // Draw field map
  // ViewField fieldView;
  // fieldView.SetComponent(&nebem);
  // fieldView.SetArea(-0.25, -1.5 * dy, -0.25, 2 * halfLX + 0.25, 1.5 * dy, 2 * halfLY + 0.25);
  // fieldView.SetPlaneXY();
  // fieldView.PlotContour("e");

  // Draw equipotential map
  ViewField potView;
  potView.SetComponent(&nebem);
  potView.SetArea(-0.25, -2. * dy, -0.25, 2 * halfLX + 0.25, 2. * dy, 2 * halfLY + 0.25);
  potView.SetPlaneXY();
  potView.PlotContour();

  Sensor sensor;
  sensor.AddComponent(&nebem);
  for (int i = 0; i < NWireX; i++)
  {
    sensor.AddElectrode(&nebem, Form("X%d", i));
  }
  for (int i = 0; i < NWireY; i++)
  {
    sensor.AddElectrode(&nebem, Form("Y%d", i));
  }
  const double tstep = 1;
  const double tmin = -0.5 * tstep;
  const unsigned int nbins = 1000;
  sensor.SetTimeWindow(tmin, tstep, nbins);

  TrackSrim track;
  track.SetSensor(&sensor);
  const std::string file = "Xenon_in_isobutane.txt";
  if (!track.ReadFile(file))
  {
    std::cerr << "Reading SRIM file failed.\n";
    return 1;
  }
  track.SetKineticEnergy(620.e6); // 4.56MeV/u
  track.SetWorkFunction(25.2);
  track.SetFanoFactor(0.25);
  track.SetAtomicMassNumbers(136, 54);
  track.SetTargetClusterSize(500);
  track.EnableTransverseStraggling(false);
  // track.SetClustersMaximum(1000);

  // DriftLineRKF drift;
  // drift.SetSensor(&sensor);
  // // Switch on signal calculation.
  // drift.EnableSignalCalculation();

  // DriftLineRKF drift;
  // drift.SetSensor(&sensor);
  // drift.EnableSignalCalculation();
  // drift.EnableAvalanche();
  // drift.EnableIonTail();

  AvalancheMicroscopic aval;
  aval.SetSensor(&sensor);
  // aval.EnableSignalCalculation();
  AvalancheMC driftMC;
  driftMC.SetSensor(&sensor);
  driftMC.SetDistanceSteps(1.e-4);
  driftMC.EnableSignalCalculation();

  TCanvas *cD = new TCanvas("cD", "cD", 600, 600);
  ViewDrift driftView;
  driftView.SetCanvas(cD);
  track.EnablePlotting(&driftView);
  // drift.EnablePlotting(&driftView);
  aval.EnablePlotting(&driftView);
  driftMC.EnablePlotting(&driftView);

  const double x0 = halfLX + 5. * rWire; // 0.25 * periodY; //2 * rWire;
  const double y0 = 1. * dy;
  const double z0 = halfLY;
  track.NewTrack(x0, y0, z0, 0, 0, -1, 0);
  // Loop over the clusters along the track.
  double xc = 0., yc = 0., zc = 0., tc = 0., ec = 0., extra = 0.;
  int nc = 0;
  int count = 0;
  while (track.GetCluster(xc, yc, zc, tc, nc, ec, extra))
  {
    count++;
    cout << "nc : " << nc << endl;
    // drift.SetElectronSignalScalingFactor(nc);
    // drift.DriftElectron(xc, yc, zc, tc);
    // drift.SetIonSignalScalingFactor(nc);
    // drift.DriftIon(xc, yc, zc, tc);
    // double ne, ni;
    // drift.GetAvalancheSize(ne, ni);
    // cout << "ne : " << ne << endl;
    // cout << "ni : " << ni << endl;

    // double xe = 0., ye = 0., ze = 0., te = 0.;
    // int ste = 0;
    // drift.GetEndPoint(xe, ye, ze, te, ste);
    // // aval.AvalancheElectron(xe + 2. * rWire, ye + 2. * rWire, ze + 2. * rWire, te, 0.1, 0, 0, 0); // 0.1 eV, random direction

    aval.AvalancheElectron(xc, yc, zc, tc, 0.1, 0, 0, 0); // 0.1 eV, random direction
    int ne = 0, ni = 0;
    // Get the number of electrons and ions in the avalanche.
    aval.GetAvalancheSize(ne, ni);
    cout << "ne : " << ne << endl;
    cout << "ni : " << ni << endl;
    // Get the number of electron drift lines.
    int np = aval.GetNumberOfElectronEndpoints();
    cout << "np : " << np << endl;
    // Initial position and time
    double x1, y1, z1, t1;
    // Final position and time
    double x2, y2, z2, t2;
    // Initial and final energy
    double e1, e2;
    // Flag indicating why the tracking of an electron was stopped.
    int status;
    // Loop over the avalanche electrons.
    for (int i = 0; i < np; ++i)
    {
      aval.GetElectronEndpoint(i, x1, y1, z1, t1, e1,
                               x2, y2, z2, t2, e2, status);
      driftMC.SetIonSignalScalingFactor(nc);
      driftMC.DriftIon(x1, y1, z1, t1);
    }

    // double xi1, yi1, zi1, ti1;
    // double xi2, yi2, zi2, ti2;
    // driftMC.GetIonEndpoint(0, xi1, yi1, zi1, ti1,
    //                      xi2, yi2, zi2, ti2, status);
    // cout << "ion t1 : " << ti1 << endl;
    // cout << "ion t2 : " << ti2 << endl;

    // Loop over the electrons in the cluster.
    // for (int k = 0; k < nc; ++k)
    // {
    //   // double xe = 0., ye = 0., ze = 0., te = 0., ee = 0.;
    //   // double dxe = 0., dye = 0., dze = 0.;
    //   // track.GetElectron(k, xe, ye, ze, te, ee, dxe, dye, dze);

    //   // double xci = 0., yci = 0., zci = 0., tci = 0., eci = 0.;
    //   // double dxe = 0., dye = 0., dze = 0.;
    //   // track.GetIon(k, xci, yci, zci, tci);
    //   // drift.DriftIon(xci, yci, zci, tci);

    //   drift.DriftElectron(xc, yc, zc, tc);
    //   drift.DriftIon(xc, yc, zc, tc);

    //   // In case of a null vector, the initial direction is randomized.
    //   // double dx0 = 0., dy0 = 0., dz0 = 0.;
    //   // Calculate an electron avalanche.
    //   // aval.AvalancheElectron(xe, ye, ze, te, ee, dxe, dye, dze);

    //   // int ne, ni;
    //   // // Get the number of electrons and ions in the avalanche.
    //   // aval.GetAvalancheSize(ne, ni);
    //   // cout << "ne : " << ne << endl;
    //   // cout << "ni : " << ni << endl;
    //   // // Get the number of electron drift lines.
    //   // int np = aval.GetNumberOfElectronEndpoints();
    //   // cout << "np : " << np << endl;
    //   // // Initial position and time
    //   // double x1, y1, z1, t1;
    //   // // Final position and time
    //   // double x2, y2, z2, t2;
    //   // // Initial and final energy
    //   // double e1, e2;
    //   // // Flag indicating why the tracking of an electron was stopped.
    //   // int status;
    //   // // Loop over the avalanche electrons.
    //   // for (int i = 0; i < np; ++i)
    //   // {
    //   //   aval.GetElectronEndpoint(i, x1, y1, z1, t1, e1,
    //   //                            x2, y2, z2, t2, e2, status);
    //   //   drift.DriftIon(x1, y1, z1, t1);
    //   // }

    //   // double xi1, yi1, zi1, ti1;
    //   // double xi2, yi2, zi2, ti2;
    //   // drift.GetIonEndpoint(0, xi1, yi1, zi1, ti1,
    //   //                      xi2, yi2, zi2, ti2, status);
    //   // cout << "ion t1 : " << ti1 << endl;
    //   // cout << "ion t2 : " << ti2 << endl;

    //   // double xe = 0., ye = 0., ze = 0., te = 0., ee = 0.;
    //   // double dx = 0., dy = 0., dz = 0.;
    //   // track.GetElectron(k, xe, ye, ze, te, ee, dx, dy, dz);
    //   // drift.DriftElectron(xe, ye, ze, te);
    // }
  }
  cout << "count : " << count << endl;

  // ViewCell cellView;
  // cellView.SetCanvas(cD);
  // cellView.SetComponent(&nebem);

  // cellView.Plot2d();
  constexpr bool twod = true;
  constexpr bool drawaxis = true;
  driftView.Plot(twod, drawaxis);

  ViewSignal signalView;
  signalView.SetSensor(&sensor);
  signalView.PlotSignal("Y2");
  ViewSignal signalView2;
  signalView2.SetSensor(&sensor);
  signalView2.PlotSignal("Y3");
  ViewSignal signalView3;
  signalView3.SetSensor(&sensor);
  signalView3.PlotSignal("Y4");
  ViewSignal signalView4;
  signalView4.SetSensor(&sensor);
  signalView4.PlotSignal("Y5");

  // std::ofstream outfile;
  // outfile.open("signal.txt", std::ios::out);
  // for (unsigned int i = 0; i < nbins; ++i)
  // {
  //   const double t = (i + 0.5) * tstep;
  //   const double f = sensor.GetSignal("Y4", i);
  //   const double fe = sensor.GetElectronSignal("Y4", i);
  //   const double fh = sensor.GetIonSignal("Y4", i);
  //   outfile << t << " " << f << " " << fe << " " << fh << "\n";
  // }

  app.Run(true);
}