#include <iostream>
#include <TApplication.h>
#include <TH1F.h>
#include <TCanvas.h>

#include "Garfield/MediumSilicon.hh"
#include "Garfield/ComponentConstant.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/TrackTrim.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/ComponentAnsys123.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/TrackHeed.hh"

using namespace Garfield;

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

  // Application
  TApplication app("app", &argc, argv);
 
  ComponentAnsys123 fm;
  fm.Initialise("ELIST.lis", "NLIST.lis", "MPLIST.lis", "PRNSOL.lis", "m");
  fm.PrintRange();
  fm.PrintMaterials();
  fm.SetWeightingField("weight.lis", "readout");
  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(-1, 0, 0, 0.5 * (x0 + x1), 0, 0);
    // fieldView.PlotContour();
    fieldView.PlotContourWeightingField("readout", "v");
  }

  // Define the medium.
  MediumMagboltz gas;
  gas.LoadGasFile("He-740Torr.gas");
  fm.SetMedium(0, &gas);

  Sensor sensor;
  sensor.AddComponent(&fm);
  sensor.AddElectrode(&fm, "readout");
  // sensor.SetTimeWindow(0., 1000., 100);
  sensor.SetTimeWindow(0., 1., 500);

  // Read the TRIM output file. 
  TrackTrim tr;
  const std::string filename = "EXYZ-He.txt";
  // Import the first 100 ions.
  if (!tr.ReadFile(filename, 500)) {
    std::cerr << "Reading TRIM EXYZ file failed.\n";
    return 1;
  }
  tr.SetWorkFunction(42);
  tr.SetFanoFactor(0.17);
  // Connect the track to a sensor.
  tr.SetSensor(&sensor);

  DriftLineRKF drift;
  drift.SetSensor(&sensor);
  drift.SetMaximumStepSize(0.01);

  // Plot the track and the drift lines.
  ViewDrift driftView;
  tr.EnablePlotting(&driftView);
  drift.EnablePlotting(&driftView);
  // drift.UseWeightingPotential();
 
  // Simulate an ion track.
  // tr.NewTrack(-11.7, 18.9427, -7.53, 0., 0., -1., 0.);
  tr.NewTrack(-11.7, 20.1, -7.53, 0., 0., -1., 0.);
  // Loop over the clusters.
  double xc, yc, zc, tc, ec, ekin;
  int ne = 0;
  while (tr.GetCluster(xc, yc, zc, tc, ne, ec, ekin)) {
    if (ne == 0) continue;
    // Simulate electron and ion drift lines starting 
    // from the cluster position. 
    // Scale the induced current by the number of electron/ion pairs 
    // in the cluster.
    std::cout << "yc = " << yc << "; ne = " << ne << "\n";
    drift.SetElectronSignalScalingFactor(ne);
    drift.DriftElectron(xc, yc, zc, tc);
  }
  driftView.SetArea(-16.704, 17.292, -6.704, 20.5927);
  driftView.Plot(true);

  ViewSignal signalView;
  signalView.SetSensor(&sensor);
  signalView.PlotSignal("readout", true, true, true);

  app.Run();
  return 0;
}