#include <iostream>
#include <fstream>

#include <TApplication.h>

#include "Garfield/SolidBox.hh" 
#include "Garfield/SolidHole.hh" 
#include "Garfield/GeometrySimple.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/MediumConductor.hh"
#include "Garfield/MediumPlastic.hh"
#include "Garfield/ComponentNeBem3d.hh"
#include "Garfield/ViewGeometry.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/Medium.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/Random.hh"
#include "Garfield/Plotting.hh"
#include "Garfield/ComponentAnalyticField.hh"

using namespace Garfield;

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

    TApplication app("app", &argc, argv);
    plottingEngine.SetDefaultStyle();

    MediumMagboltz gas;

    gas.LoadGasFile("ArIso5new.gas");
//    gas.LoadGasFile("ArIso5Oxygen0p1.gas");
//    gas.LoadGasFile("ArIso5_1pcO2.gas");

    const std::string path = std::getenv("GARFIELD_INSTALL");
  	gas.LoadIonMobility(path + "/share/Garfield/Data/IonMobility_Ar+_Ar.txt");
  	gas.LoadIonMobility(path + "/share/Garfield/Data/IonMobility_C8Hn+_iC4H10.txt");

    gas.EnableDrift(true);
    gas.EnableThermalMotion(true);
    gas.EnablePrimaryIonisation(true);

  // Set the Penning transfer efficiency.
  constexpr double rPenning = 0.26; // "3D simulation of micromegas detector performance" GUO Jun-Jun
  constexpr double lambdaPenning = 0.;
    gas.EnablePenningTransfer(rPenning, lambdaPenning, "ar");
    gas.SetMaxElectronEnergy(200.);
    gas.Initialise(true);
 




  const double xmin = -0.01;
  const double xmax =  0.01;
  const double ymin = -0.002;  // 0.001 ORG
  const double ymax =  0.008;


    // 50 um amplification gap.
    ComponentAnalyticField cmpA;
    cmpA.SetMedium(&gas);
    constexpr double yMesh = 0.005;
    constexpr double vMesh = 0.;
    constexpr double vPad = 320.; // 320 org
    cmpA.AddPlaneY(yMesh, vMesh, "mesh1");
    cmpA.AddPlaneY(0., vPad, "pad");

    Sensor sensorA;
    sensorA.AddComponent(&cmpA);
    sensorA.SetArea(-5.0, 0., -5.0, 5.0, yMesh, 5.0);

// avalanche: field, eta
std::cout << "\nField intensity in the middle of avalanche region:\n";
double y2 = 0.5 * (yMesh); 
auto edrift2 = cmpA.ElectricField(0., y2, 0.);
double magnitude2 = sqrt(edrift2[0] * edrift2[0] + edrift2[1] * edrift2[1]);
std::cout << "At y2 = " << y2 << " cm:\n";
std::cout << "  E = " << magnitude2/1000 << " kV/cm\n";
double eta2 = 0.;
gas.ElectronAttachment(edrift2[0], edrift2[1], edrift2[2], 0., 0., 0., eta2);
std::cout <<" eta2 = " << eta2 << "\n";

  ViewField fieldView;
  
  fieldView.SetComponent(&cmpA);
  fieldView.SetPlaneXY();
  // Set the plot limits in the current viewing plane.

  fieldView.SetArea(xmin, ymin, xmax, ymax);
  fieldView.SetVoltageRange(0., 400.);
  fieldView.SetNumberOfContours(10);
  fieldView.PlotContour();
    std::vector<double> xf;
    std::vector<double> yf;
    std::vector<double> zf;

 fieldView.EqualFluxIntervals(xmin, 0.0048, 0, xmax, 0.0048, 0, xf, yf, zf, 100); // 0.0056,
 fieldView.PlotFieldLines(xf, yf, zf, true, false); 



    AvalancheMC avalMC;
    avalMC.SetSensor(&sensorA);

   avalMC.SetDistanceSteps(1.e-4);
//   avalMC.EnableAvalancheSizeLimit(10);
//   avalMC.SetCollisionSteps(100);
//   avalMC.UseWeightingPotential(false);
//   avalMC.DisableAvalancheSizeLimit();


  ViewDrift driftView;
  constexpr bool plotDrift = true;  // true or false  SOS
  if (plotDrift) {
    avalMC.EnablePlotting(&driftView);
  }

  unsigned int nTotal = 0;
  unsigned int nBF = 0;
  constexpr unsigned int nEvents = 1;
  for (unsigned int i = 0; i < nEvents; ++i) { 

      double x0 = 0.0; 
      double y0 = yMesh - 0.0001; //0.0045; 
      double z0 = 0.0; //0.005 - 1e-14;
      double t0 = 0.;
      double e0 = 0.1;
        // The above is with the initial electron travelling down (direction given by last 3 coordinates)      
        //aval.AvalancheElectron(xi, yi, zi, ti, ei, 0, 0, -1);
      double xee = 0., yee = 0., zee = 0.0045, tee = 0.;
      int status = 0;

  int ElectronNumber=0;


//avalMC.DriftElectron(x0, y0, z0, 0.0);
avalMC.AvalancheElectron(x0, y0, z0, 0.0, false);
// avalMC.AvalancheElectron(x2, yMesh - 1e-14, z2, t2, false);




const unsigned int np = avalMC.GetNumberOfElectronEndpoints();
    std::cout << np << "/" << nEvents << "\n";

  }
 

  if (plotDrift) {
    TCanvas* c2 = new TCanvas();
 {
      driftView.SetPlane(0, 0, -1, 0, 0, 0);
      driftView.SetArea(xmin, ymin, xmax, ymax); //(-2 * wirepitch, -0.02, 2 * wirepitch, 0.02);
      driftView.SetCanvas(c2);
      constexpr bool twod = true;
      driftView.Plot(twod);
    }
  }

/////////////////////////////////////////////////////////////////////////////


  
 app.Run(true);
}