#include <cstdlib>
#include <iostream>
#include <fstream>

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

#include "Garfield/ComponentAnsys123.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ViewFEMesh.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/Random.hh"
#include "Garfield/SolidBox.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/Plotting.hh"
#include "Garfield/ViewGeometry.hh"
#include "Garfield/ViewCell.hh"
#include "Garfield/DriftLineRKF.hh"


using namespace Garfield;
using namespace std;

int main(int argc, char * argv[]) {
  TApplication app("app", &argc, argv);

  // Setup the gas.
  MediumMagboltz gas("ar", 70., "co2", 29., "o2", 1.0);
  gas.SetTemperature(293.15);
  gas.SetPressure(760.);
  gas.Initialise(true);
  
  // Set the Penning transfer efficiency.
  constexpr double rPenning = 0.51;
  constexpr double lambdaPenning = 0.;
  gas.EnablePenningTransfer(rPenning, lambdaPenning, "ar");
  // Load the ion mobilities.
  const std::string path = std::getenv("GARFIELD_INSTALL");
  gas.LoadIonMobility(path + "/share/Garfield/Data/IonMobility_Ar+_Ar.txt");

  // Load the field map.
  ComponentAnsys123 fm;
  fm.Initialise("ELIST.lis", "NLIST.lis", "MPLIST.lis", "PRNSOL.lis", "mm");
  fm.EnableMirrorPeriodicityX();
  fm.EnableMirrorPeriodicityY();
  fm.PrintRange();

  // Associate the gas with the corresponding field map material.
  fm.SetGas(&gas); 
  fm.PrintMaterials();
  //fm.Check();

  // Dimensions of the GEM [cm]
  constexpr double pitch = 0.014;
  
  const double width = 2.;
  ComponentAnalyticField* comp = new ComponentAnalyticField();
  comp->AddPlaneY(-width / 2, 0., "Cathode");
  comp->AddPlaneY(width / 2, 1000., "Anode");
  comp->SetMagneticField(0.5, 0., 0);
 
  // Create the sensor.
  Sensor sensor;
  sensor.AddComponent(&fm);
  sensor.SetArea(-5.0 * pitch, -5.0 * pitch, -0.02,
                  5.0 * pitch,  5.0 * pitch,  0.025);

  AvalancheMicroscopic aval;
  aval.SetSensor(&sensor);

  AvalancheMC drift;
  drift.SetSensor(&sensor);
  drift.SetDistanceSteps(2.e-4);

  ViewDrift driftView;
  ViewField fieldView;
  ViewFEMesh meshView;
  constexpr bool plotDrift = true;
  if (plotDrift) {
    aval.EnablePlotting(&driftView);
    drift.EnablePlotting(&driftView);
  }
  int netot = 0;
  for (unsigned int j = 0; j < 1 ; ++j) {
    constexpr unsigned int nEvents = 5;
    for (unsigned int i = 0; i < nEvents ; ++i) { 
      // Randomize the initial position.
      const double x0 = 0. * pitch;
      const double y0 = -0.5 * pitch + RndmUniform() * pitch;
      const double z0 = 0.02; 
      const double t0 = 0.;
      const double e0 = 0.1;
      aval.AvalancheElectron(x0, y0, z0, t0, e0, 0., 0., 0.);
      int ne = 0, ni = 0;
      aval.GetAvalancheSize(ne, ni);
      netot += ne;
      const unsigned int np = aval.GetNumberOfElectronEndpoints();
      double xe1, ye1, ze1, te1, e1;
      double xe2, ye2, ze2, te2, e2;
      int status;
      for (unsigned int k = 0; k < np; ++k) {
        aval.GetElectronEndpoint(k, xe1, ye1, ze1, te1, e1, 
                                    xe2, ye2, ze2, te2, e2, status);
      }
    cout << "xe2:=" << xe2 << endl;
    cout << "ye2:=" << ye2 << endl; 
    cout << "ze2:=" << ze2 << endl; 
    cout << "ne: " << ne << endl;       
    }
  cout << "total electron: " << netot << endl;
  cout << "total gain:= " << netot / 5.000 << endl;
  }
  if (plotDrift) {
    TCanvas* cd = new TCanvas();
    constexpr bool plotMesh = true;
    if (plotMesh) {
      meshView.SetCanvas(cd);
      meshView.SetComponent(&fm);
      constexpr bool twod = true;
      // x-z projection.
      meshView.SetPlane(0, -1, 0, 0, 0, 0);
      if (twod) {
        meshView.SetArea(-2 * pitch, -0.02, 2 * pitch, 0.02);
      } else {
        meshView.SetArea(-0.5 * pitch, -0.5 * pitch, -0.02, 0.5 * pitch, 0.5 * pitch, 0.02);
      }
      meshView.SetFillMesh(true);
      meshView.SetColor(0, kGray);
      // Set the color of the kapton.
      meshView.SetColor(2, kYellow + 3);
      meshView.EnableAxes();
      meshView.SetViewDrift(&driftView);
      meshView.Plot(twod);
    } else {
      driftView.SetPlane(0, -1, 0, 0, 0, 0);
      driftView.SetArea(-2 * pitch, -0.02, 2 * pitch, 0.02);
      driftView.SetCanvas(cd);
      constexpr bool twod = true;
      driftView.Plot(twod);
    }
  }
  cout << "Done2.\n";
  app.Run();
}