#include <iostream>
#include <TROOT.h>
#include <TApplication.h>
#include <math.h>
#include <bits/stdc++.h>
#include <omp.h>
#include <vector>
#include <fstream>
#include <cstddef> // for std::size_t
#include <stdio.h> // for log files

#include "Garfield/MediumMagboltz.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/FundamentalConstants.hh"
#include "Garfield/Random.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/ViewCell.hh"
#include "Garfield/ViewDrift.hh"
using namespace Garfield;

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

  // https://root-forum.cern.ch/t/rewriting-macro-as-c-application/3225
  int argcROOT = 1;
  char* argvROOT = argv[0];
  TApplication theApp("App", &argcROOT, &argvROOT);

  time_t tstart,end;
  double time_taken;
  time(&tstart);
  
  // Wire radii [cm].
  // Anode Wire
  const double rWire = .008;
  // Guard wire
  const double rTube = 1.0;

  // Setup the gas.
  MediumMagboltz gas;
  // Set the temperature [K] and pressure [Torr].
  gas.SetTemperature(293.15);
  gas.SetPressure(50.);
  // Set the gas mixture.
  gas.SetComposition("isobutane", 100.);
  // Read the ion mobility table from file.
  gas.LoadIonMobility("isobutane_ion.txt");
  gas.SetMaxElectronEnergy(3000.);
  gas.Initialise(true);
  
  // Setup the electric field.
  ComponentAnalyticField cmp;
  cmp.SetMedium(&gas);
  
  // Add anode wire
  cmp.AddWire(0., 0., 2.*rWire, 1000., "anode");
  // Add outer cylinder
  cmp.AddTube(rTube, 0., 0., "tube");

  // Weighting field activation
  cmp.AddReadout("anode");

  // Make a sensor.
  Sensor sensor;
  sensor.AddComponent(&cmp);
  sensor.AddElectrode(&cmp, "anode");
  
  // Set the time bins.
  const unsigned int nTimeBins = 1000;
  const double tmin =  0.;
  const double tmax = 1000.;
  const double tstep = (tmax - tmin) / nTimeBins;
  sensor.SetTimeWindow(tmin, tstep, nTimeBins);

  TrackHeed track;
  track.SetSensor(&sensor);
  track.SetParticle("proton");
  track.SetKineticEnergy(10.e6);
  
  double xc = 0., yc = 0., zc = 0., tc = 0., ec = 0., extra = 0.;
  int nc = 0;
  std::vector<double> electrons;
  double xe = 0., ye = 0., ze = 0., te = 0., ee = 0.;
  double dx = 0., dy = 0., dz = 0.;
  int NofProton = 100;
  double rTrack = 0.5;
  double x0 = rTrack;
  double y0 = -sqrt(rTube * rTube - rTrack * rTrack);
  int index = 0;
  for (int i = 0; i < NofProton; i++) {
    // Simulate a track
    track.NewTrack(x0, y0, 0, 0, 0, 1, 0);
		while (track.GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
			for (int k = 0; k < nc; ++k) {
			  track.GetElectron(k, xe, ye, ze, te, ee, dx, dy, dz);
			  electrons.push_back(xe);
			  electrons.push_back(ye);
			  electrons.push_back(ze);
			  electrons.push_back(te);
			  ++index;
			}
		}
  }

  double xe1, ye1, ze1, te1, e1;
  double xe2, ye2, ze2, te2, e2;
  int status;
  bool calculate_signal = true;
  AvalancheMicroscopic aval;
  aval.SetSensor(&sensor);
  aval.EnableSignalCalculation(calculate_signal);
  DriftLineRKF drift;
  drift.SetSensor(&sensor);
  drift.EnableSignalCalculation(calculate_signal);
  constexpr bool plotDrift = false;
  TCanvas* cD = nullptr;
  ViewCell cellView;
  ViewDrift driftView;
  if (plotDrift) {  	
    cD = new TCanvas("cD", "", 600, 600);
    cellView.SetCanvas(cD);
    cellView.SetComponent(&cmp);
    driftView.SetCanvas(cD);
    aval.EnablePlotting(&driftView);
    drift.EnablePlotting(&driftView);
  }
  std::cout << "Total electrons: " << index << std::endl;
  for (int k = 0; k < index; k++) { 
    aval.AvalancheElectron(electrons.at(k*4 + 0),electrons.at(k*4 + 1),electrons.at(k*4 + 2),electrons.at(k*4 + 3), 0.1, 0,0,0);
    drift.DriftIon(electrons.at(k*4 + 0),electrons.at(k*4 + 1),electrons.at(k*4 + 2),electrons.at(k*4 + 3));
    const int np = aval.GetNumberOfElectronEndpoints();
    std::cout << k << std::endl;
    for (int j = np; j--;) {
      aval.GetElectronEndpoint(j, xe1, ye1, ze1, te1, e1, 
                                  xe2, ye2, ze2, te2, e2, status);
      drift.DriftIon(xe1, ye1, ze1, te1);
    }
  }

	// Plot the cell layout and drift lines if requested
	if (plotDrift) {
    cD->Clear();
    cellView.SetComponent(&cmp);
	cellView.SetArea(-2., -2., 2., 2.);
    cellView.Plot2d();
    constexpr bool twod = true;
    constexpr bool drawaxis = false;
    driftView.Plot(twod, drawaxis);
  }

  // Write things into files
  double t = 0;

  // Write raw signals of non-pad electrodes
  std::ofstream total;
  total.open("tot_signals.txt", std::ios::out);
  std::ofstream ion;
  ion.open("ion_signals.txt", std::ios::out);
  std::ofstream el;
  el.open("el_signals.txt", std::ios::out);
  //first line for titles of columns
  total << "Time" << "\t" << "Anode" <<  std::endl;
  ion << "Time" << "\t" << "Anode" <<  std::endl;
  el << "Time" << "\t" << "Anode" <<  std::endl;
  for (unsigned int j = 0; j < nTimeBins; ++j) {
		t = (j + 0.5) * tstep;
		total << t << "\t" << sensor.GetSignal("anode",j) << std::endl;
		ion << t << "\t" << sensor.GetIonSignal("anode",j) << std::endl;
		el << t << "\t" << sensor.GetElectronSignal("anode",j) << std::endl;
  }
  total.close();
  ion.close();
  el.close();
	
  time(&end);
  time_taken = double(end-tstart);
  std::cout << "Time taken by program is: " << std::fixed << time_taken << std::setprecision(5) << " sec" << std::endl;
  if (plotDrift) theApp.Run(true);
}