#include <iostream>
#include <fstream>
#include <cstdlib>

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

#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/TrackSrim.hh"
#include "Garfield/ViewDrift.hh"

using namespace Garfield;

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

  TApplication app("app", &argc, argv);
 
  // Make a gas medium.
  MediumMagboltz gas;
  gas.LoadGasFile("Ar90_CH4_10_1bar.gas");
  gas.LoadIonMobility("IonMobility_Ar+_Ar.txt");

  // Make a component with analytic electric field.
  ComponentAnalyticField cmp;
  cmp.SetMedium(&gas);
  // Wire radius [cm]
  const double rWire = 25.e-4;
  // Outer radius of the tube [cm]
  const double rTube = 0.71;
  // Voltages
  const double vWire = 2730.;
  const double vTube = 0.;
  // Add the wire in the centre.
  cmp.AddWire(0, 0, 2 * rWire, vWire, "s");
  // Add the tube.
  cmp.AddTube(rTube, vTube, 0);

  // Make a sensor.
  Sensor sensor(&cmp);
  sensor.AddElectrode(&cmp, "s");

  ViewDrift driftView;

  TrackSrim track(&sensor);
  track.ReadFile("Hydrogen_in_Ar_90_CH4_10.txt");
  track.SetKineticEnergy(2.1e6);  

  AvalancheMC drift(&sensor);
  drift.SetDistanceSteps(0.01);
 
  TCanvas* cD = nullptr;
  constexpr bool plotDrift = false; 
  if (plotDrift) {
    cD = new TCanvas("cD", "", 600, 600);
    driftView.SetCanvas(cD);
    drift.EnablePlotting(&driftView);
    track.EnablePlotting(&driftView);
  }

  const double rTrack = 0.3;
  const double x0 = rTrack;
  const double y0 = -sqrt(rTube * rTube - rTrack * rTrack);
  const unsigned int nTracks = 100;
  for (unsigned int j = 0; j < nTracks; ++j) {
    track.NewTrack(x0, y0, 0, 0, 0, 1, 0);
    std::cout << "There are " << track.GetClusters().size() 
              << " clusters on the track.\n";
    for (const auto& cluster : track.GetClusters()) {
      const int num_electrons = cluster.n;
      std::cout << "  Cluster with " << num_electrons << " electrons.\n";
      for (int jj = 0; jj < num_electrons; jj++) {
        drift.DriftElectron(cluster.x, cluster.y, cluster.z, cluster.t);
        drift.DriftIon(cluster.x, cluster.y, cluster.z, cluster.t);
      }
    }
  }
  std::cout << "Done.\n";
  if (plotDrift) {
    cD->Clear();
    cmp.PlotCell(cD);
    constexpr bool twod = true;
    constexpr bool drawaxis = false;
    driftView.Plot(twod, drawaxis);
  }
  app.Run();
}