#include <iostream>
#include <fstream>

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

#include "Garfield/AvalancheMC.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Plotting.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/TrackSrim.hh"
#include "Garfield/ViewSignal.hh"

using namespace Garfield;

int main(int argc, char * argv[]) {
  TApplication app("app", &argc, argv);
  
  double ratioofAr = 0.85;
  double ratioofCO2 = 1. - ratioofAr;

  MediumMagboltz gas;
  const double pressure = 0.26;
  const double temperature = 300;
  gas.SetComposition("Ar", ratioofAr, "CO2", ratioofCO2);
  gas.SetTemperature(temperature);
  gas.SetPressure(pressure*AtmosphericPressure);                        
  
  gas.SetMaxElectronEnergy(600);
  gas.Initialise(true);  
  gas.LoadIonMobility("IonMobility_Ar+_Ar.txt");

  ComponentAnalyticField cmp;
  cmp.SetMedium(&gas);
  const double rWire = 15.e-4;
  const double rTube = 3.5;
  const double vWire = 1100;
  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);

  Sensor sensor(&cmp);
  sensor.AddElectrode(&cmp, "s");
  // Set the signal time window.
  const double tstep = 0.2;
  const double tmin = 0;
  const double tmax = 10000.;
  const unsigned int nbins = (tmax - tmin) / tstep;
  sensor.SetTimeWindow(tmin, tstep, nbins);

  TrackSrim tr(&sensor);
  const std::string file = "Li_in_Ar_CO2_gas.txt";
  if (!tr.ReadFile(file)) {
    std::cerr << "Reading SRIM file failed.\n";
  }
  tr.SetKineticEnergy(0.84e6);
  // tr.SetTargetClusterSize(10);
  tr.EnableTransverseStraggling(false);
  tr.EnableLongitudinalStraggling(false);
  tr.Print();

  const double rTrack = 0.3;
  const double projectiley = -sqrt(rTube * rTube - rTrack * rTrack);
  if (!tr.NewTrack(rTrack, projectiley, 0., 0., 0., 1., 0.)) {
    std::cerr << "Generating clusters failed.\n";
    return 0;
  }
  std::cout << tr.GetClusters().size() << " clusters on the track.\n";
  for (const auto& cluster : tr.GetClusters()) {
    AvalancheMicroscopic driftelectron(&sensor);
    AvalancheMC driftion(&sensor);
    driftion.SetDistanceSteps(2.e-4);
    std::cout << "Cluster with " << cluster.n << " electron-ion pairs.\n";
    for (int k = 0; k < cluster.n; ++k) {
      std::cout << "  Electron " << k << ":\n";
      driftelectron.AvalancheElectron(cluster.x, cluster.y, cluster.z, 
                                      cluster.t, 0.01, 0., 0., 0.);
      int ne = 0, ni = 0;
      driftelectron.GetAvalancheSize(ne, ni);
      std::cout << "    " << ne << " electrons, " << ni << " ions.\n";
      // Drift the ions created in the avalanche.
      for (const auto& electron : driftelectron.GetElectrons()) {
        const auto& p0 = electron.path[0];
        driftion.DriftIon(p0.x, p0.y, p0.z, p0.t);       
      }
    }
  }

  TCanvas* cS = new TCanvas("cS", "", 600, 600);
  ViewSignal signalView(&sensor);
  signalView.SetCanvas(cS);
  signalView.PlotSignal("s", "tei");
  std::cout << "Calculation is finished.\n";
  app.Run();
}