#include <iostream>
#include <cstdlib>

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

#include "Garfield/FundamentalConstants.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/ComponentAnsys123.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ViewFEMesh.hh"
#include "Garfield/ViewMedium.hh"

using namespace Garfield;

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

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

  ComponentAnsys123 fm;
  fm.Initialise("ELIST.lis", "NLIST.lis", "MPLIST.lis", "PRNSOL.lis", "cm");
  fm.PrintRange();
  fm.PrintMaterials();
  // Get the extent of the field map.
  double x0, y0, z0, x1, y1, z1;
  fm.GetBoundingBox(x0, y0, z0, x1, y1, z1);

  ViewField fieldView;
  constexpr bool plotField = true;
  if (plotField) {
    fieldView.SetComponent(&fm);
    // Set the normal vector of the viewing plane.
    fieldView.SetPlane(0, 0, -1, 0, 0, 0.5 * (z0 + z1));
    fieldView.PlotContour();
  }

  Sensor sensor;
  sensor.AddComponent(&fm);

  MediumMagboltz gas;
  // gas.LoadGasFile("4He.gas");
  gas.LoadGasFile("ar_93_co2_7.gas");

  fm.SetMedium(0, &gas);
  fm.PrintMaterials();

  TrackHeed track;
  track.SetParticle("proton");
  track.SetEnergy(100.e9);
  track.SetSensor(&sensor);
  track.Initialise(&gas, true);

  DriftLineRKF drift;
  drift.SetSensor(&sensor);
  
  ViewDrift driftview;
  TCanvas* CD = nullptr;
  constexpr bool plotDrift = true;
  if (plotDrift) {
    CD = new TCanvas("CD", "", 600, 600);
    driftview.SetCanvas(CD);
    drift.EnablePlotting(&driftview);
    track.EnablePlotting(&driftview);
  }

  const unsigned int nTracks = 1;
  size_t nClusters = 0;
  for (unsigned int j = 0; j < nTracks; ++j) {
    sensor.ClearSignal();
    // Starting point of the track.
    double xt = -0.1, yt = 0.204, zt = -0.07, tt = 0.;
    // Initial direction of the track.
    double dxt = 0., dyt = -1., dzt = 0.; 
    track.NewTrack(xt, yt, zt, tt, dxt, dyt, dzt);
    // Loop over the clusters on the track.
    double xc = 0., yc = 0., zc = 0., tc = 0., ec = 0., extra = 0.;
    int nc = 0;
    while (track.GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
      ++nClusters;
      // 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);
        drift.DriftElectron(xe, ye, ze, te);
      }
    }
  }
  std::cout << nClusters << " clusters...\n";
  const bool twod = true;
  const bool axis = true;
  driftview.Plot(twod, axis);
  app.Run(true);
}