// multiple tracks and calculate the average collected charge
#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-560newTorr.gas");
  fm.SetMedium(0, &gas);
  // Load the ion mobility.
  //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_He+_He.txt");


  Sensor sensor;
  sensor.AddComponent(&fm);
  sensor.AddElectrode(&fm, "readout");
  // Set the time bins for the induced current.
  const unsigned int nTimeBins = 500;
  const double tmin =  0.;
  //onst double tmax = 2000.;
  const double tmax = 500.;
  const double tstep = (tmax - tmin) / nTimeBins;
  sensor.SetTimeWindow(tmin, tstep, nTimeBins);


  // Read the TRIM output file. 
  TrackTrim tr;
  const std::string filename = "EXYZ-560He.txt";
  // Import the first 100 ions.
  if (!tr.ReadFile(filename)) {
    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();
 
  constexpr unsigned int nTracks = 1000;
  unsigned int nesum = 0;
  for (unsigned int i = 0; i < nTracks; ++i) {
    sensor.NewSignal();
    // Simulate an ion track.
    tr.NewTrack(-11.7, 20.4, -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)) {
      // Sum up the number of electron/ion pairs created.
      nesum += ne;
      // 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.
      drift.SetElectronSignalScalingFactor(ne);
      drift.DriftElectron(xc, yc, zc, tc);
      //drift.SetIonSignalScalingFactor(ne);
      //drift.DriftIon(xc, yc, zc, tc);
    }
  }
  std::cout << "Average number of electron/ion pairs: " 
            << double(nesum) / nTracks << "\n";
  sensor.IntegrateSignals();
  double qsum = sensor.GetSignal("readout", nTimeBins - 1);
  std::cout << "Integrated charge: " << qsum / ElementaryCharge << " e-\n"; 
  ViewSignal signalView;
  signalView.SetSensor(&sensor);
  // signalView.PlotSignal("readout", true, true, true);
  signalView.PlotSignal("readout", true, true);
  app.Run();
  return 0;
}