#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/DriftLineRKF.hh"
#include "Garfield/AvalancheMicroscopic.hh"

using namespace Garfield;

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

  // TApplication app("app", &argc, argv);
  // Make a gas medium.
  MediumMagboltz gas;
  gas.LoadGasFile("isobutane_50_Torr.gas");
  const std::string path = std::getenv("GARFIELD_INSTALL");
  gas.LoadIonMobility(path + "/share/Garfield/Data/IonMobility_C8Hn+_iC4H10.txt");

  // Describe the cell layout.
  ComponentAnalyticField cmp;
  cmp.SetMedium(&gas);
  // Wire radius [cm]
  const double rWire = 80.e-4;
  // Outer radius of the tube [cm]
  const double rTube = 1.;
  // Voltages
  const double vWire = 1000.;
  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, "t");
  // Request calculation of the weighting field. 
  cmp.AddReadout("s");

  // Make a sensor.
  Sensor sensor;
  sensor.AddComponent(&cmp);
  sensor.AddElectrode(&cmp, "s");
  // Set the signal time window.
  const double tstep = 0.01;
  const double tmin = -0.5 * tstep;
  const unsigned int nTimeBins = 50000;
  sensor.SetTimeWindow(tmin, tstep, nTimeBins);

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

  // Microscopic simulation.
  AvalancheMicroscopic aval;
  aval.SetSensor(&sensor);
  aval.EnableSignalCalculation();
 
  const double rTrack = 0.5;
  const double x0 = rTrack;
  const double y0 = 0.;
  const bool microscopic = true;
  if (microscopic) {
    aval.AvalancheElectron(x0, y0, 0., 0., 0.1);
    auto np = aval.GetNumberOfElectronEndpoints();
    std::cout << np << " electrons\n";
    for (unsigned int i = 0; i < np; ++i) {
      double x1, y1, z1, t1, e1;
      double x2, y2, z2, t2, e2;
      int status = 0;
      aval.GetElectronEndpoint(i, x1, y1, z1, t1, e1,
                                  x2, y2, z2, t2, e2, status);
      drift.DriftIon(x1, y1, z1, t1);
    }
  } else {
    drift.DriftElectron(x0, y0, 0., 0.);
    double ne = 0., ni = 0.;
    drift.GetAvalancheSize(ne, ni);
    std::cout << ne << " electrons\n";
  }
  std::cout << "Done.\n";
  std::ofstream outfile;
  outfile.open("signals.txt", std::ios::out);
  outfile << "#Time  Total  Electrons  Ions\n";
  for (unsigned int j = 0; j < nTimeBins; ++j) {
    double t = (j + 0.5) * tstep;
    outfile << t << "  " 
            << sensor.GetSignal("s", j) << "  "
            << sensor.GetIonSignal("s", j) << "  "
            << sensor.GetElectronSignal("s", j) << "\n";
  }
  // app.Run(true);
}