#include <iostream>
#include <TROOT.h>
#include <TApplication.h>
#include <math.h>
#include <bits/stdc++.h>
#include <omp.h>
#include <vector>

#include "Garfield/MediumMagboltz.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ViewCell.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/FundamentalConstants.hh"
#include "Garfield/Random.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/ViewFEMesh.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/AvalancheMC.hh"


using namespace Garfield;

// Impulse response of the PASA.
double transfer(double t) {
  constexpr double tau = 300;//160;
  constexpr double fC_to_mV = .0023385;//12.7;
  return fC_to_mV * exp(4) * pow((t / tau), 4) * exp(-4 * t / tau);  
}

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

  TApplication app("app", &argc, argv);

  // Wire radii [cm].
  // Active wire
  const double rActive = .00254; //1 mil
  // Guard wire
  const double rGuard = .0127; //5 mil

  // Setup the gas.
  MediumMagboltz gas;
  // Set the temperature [K] and pressure [Torr].
  gas.SetTemperature(293.15);
  gas.SetPressure(750.);
  // Set the gas mixture.
  gas.SetComposition("ne", 85.72, "co2", 9.52, "n2", 4.76);
  // Read the electron transport coefficients from a .gas file.
  // gas.LoadGasFile("Ne_90_CO2_10_N2_5_with_mg.gas");
  // Read the ion mobility table from file.
  const std::string garfpath = std::getenv("GARFIELD_HOME");
  gas.LoadIonMobility(garfpath + "/Data/IonMobility_Ne+_Ne.txt");
  gas.SetMaxElectronEnergy(200.);
  gas.Initialise(true);
 
  // Setup the electric field.
  ComponentAnalyticField cmp;
  cmp.SetMedium(&gas);
  
  // Adding the cathode plate.
  cmp.AddPlaneY(0., -641.67, "pad_plane");

  //Adding ten rectangular wire 'loops'.
  constexpr double x = 3.81; //the depth is 3 in. (~7.62 cm) centered at 0
  cmp.AddWire(-x, 0.3, 2.*rActive, -583.34, "rl");cmp.AddWire(x, 0.3, 2.*rActive, -583.34, "rl");
  cmp.AddWire(-x, 0.554, 2.*rActive, -525.01, "rl");cmp.AddWire(x, 0.554, 2.*rActive, -525.01, "rl");
  cmp.AddWire(-x, 0.808, 2.*rActive, -466.68, "rl");cmp.AddWire(x, 0.808, 2.*rActive, -466.68, "rl");
  cmp.AddWire(-x, 1.062, 2.*rActive, -408.35, "rl");cmp.AddWire(x, 1.062, 2.*rActive, -408.35, "rl");
  cmp.AddWire(-x, 1.316, 2.*rActive, -350.02, "rl");cmp.AddWire(x, 1.316, 2.*rActive, -350.02, "rl");
  cmp.AddWire(-x, 1.570, 2.*rActive, -291.69, "rl");cmp.AddWire(x, 1.570, 2.*rActive, -291.69, "rl");
  cmp.AddWire(-x, 1.824, 2.*rActive, -233.36, "rl");cmp.AddWire(x, 1.824, 2.*rActive, -233.36, "rl");
  cmp.AddWire(-x, 2.078, 2.*rActive, -175.03, "rl");cmp.AddWire(x, 2.078, 2.*rActive, -175.03, "rl");
  cmp.AddWire(-x, 2.332, 2.*rActive, -116.70, "rl");cmp.AddWire(x, 2.332, 2.*rActive, -116.70, "rl");
  cmp.AddWire(-x, 2.586, 2.*rActive, -58.37, "rl");cmp.AddWire(x, 2.586, 2.*rActive, -58.37, "rl");
  
  //Adding Frisch grid
  double spacing = .0725714286;
  int NofWires = int(7.62/spacing);
  for (int i = 0; i <= NofWires; i++) {
  cmp.AddWire(i*spacing-x, 2.84, 2.*rActive, 0., "f");
  }
  
  //Adding the anode wires and the guard wires
  //Left side
  constexpr double vG = 1000.;
  constexpr double vA = 1400.;
  const double center_left = -x + 0.5*1.5*2.54;
  const double center_right = x - 0.5*1.5*2.54;
  cmp.AddWire(center_left - .8, 3.793, 2.*rActive, vA, "al-1");
  cmp.AddWire(center_left     , 3.793, 2.*rActive, vA, "al");
  cmp.AddWire(center_left + .8, 3.793, 2.*rActive, vA, "al1");
  //Right side
  cmp.AddWire(center_right - .8, 3.793, 2.*rActive, vA, "ar-1");
  cmp.AddWire(center_right     , 3.793, 2.*rActive, vA, "ar");
  cmp.AddWire(center_right + .8, 3.793, 2.*rActive, vA, "ar+1");
  
  //"Square" blocks on the sides of anode wires
  const double dBlock = .25*2.54;
  cmp.AddWire( -x + 0.5*dBlock, 3.793, dBlock, 0., "Block");
  cmp.AddWire(-0.5*dBlock - .1, 3.793, dBlock, 0., "Block");
  cmp.AddWire( 0.5*dBlock + .1, 3.793, dBlock, 0., "Block");
  cmp.AddWire(  x - 0.5*dBlock, 3.793, dBlock, 0., "Block");
  
  
  //Adding neutral plane above the wires
  cmp.AddPlaneY(4.211, 0., "top_plane");

  // Make a sensor.
  Sensor sensor;
  sensor.SetArea(-x, 0, -40., x, 4.3, 40.); //setting drift region
  sensor.AddComponent(&cmp);
  sensor.SetTimeWindow(0., 1, 10000); 
  
  cmp.AddReadout("ar");
  sensor.AddElectrode(&cmp, "ar");
  
  TrackHeed track;
  track.SetSensor(&sensor);
  track.SetParticle("proton");
  track.SetKineticEnergy(100.e12);
  //track.SetKineticEnergy(10.e6);
  //Simulate a track
  track.NewTrack(x, 1.0, 0., 100., -1., 0., 0.);

  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 index = 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;
  double xi1, yi1, zi1, ti1;
  double xi2, yi2, zi2, ti2;
  int status;
  int size = electrons.size()/4;
  if (size != index) {
  	std::cout << "Something is wrong with electrons!!" << std::endl;
  	return 0;
  }
  #pragma omp parallel for
  for (int k = 0; k < index; k++) {
    AvalancheMicroscopic aval;
    aval.SetSensor(&sensor);
    aval.EnableSignalCalculation();
 
    aval.AvalancheElectron(center_right, 1.0, 0.0, 100., 0.1, 0,0,0);
    const int np = aval.GetNumberOfElectronEndpoints();
    DriftLineRKF drift;
    drift.SetSensor(&sensor);
    drift.EnableSignalCalculation();
    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);
    }
  }
  
  app.Run(true);
}