#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"

using namespace Garfield;

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

//  Define gas medium

    MediumMagboltz gas;

    gas.LoadGasFile("ArIso5new.gas"); // produced by magboltz of garfield++
 
    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);
 
 
  MediumConductor Cu;
  MediumPlastic Kp;
  Kp.SetDielectricConstant(4.0);

  // Define geometry
  GeometrySimple geo;
  geo.SetMedium(&gas);
  
  double drftgap = 1.0;  // cm
  double wiredia = 0.004;
  double holedia = 0.006;
  double wirepitch = wiredia + holedia;
  double andplz = 0.0;
  double andplLZ = 0.0010;
  double thickZ1 = 0.0005;
  double thickZ = 0.0045;
  double wirez = 0.00275;
  double wirez2 = 0.00525;
  double lenLX = wirepitch;
  double lenLY = wirepitch;
  double drftplLZ = 0.0010;
  double drftplz = wirez + wiredia / 2 + drftgap + drftplLZ / 2;
  double drftplV = -300.0; // -250.0
  double andplV = 320.0;  // 345.0
  double wireV1 = 0.0;

  constexpr double pitch = 0.02;
  const double xmin = -0.01;
  const double xmax =  0.01;
  const double zmin = -0.002;  // 0.001 ORG
  const double zmax =  0.008;

  // Drift plane
  SolidBox cathode(0.0, 0.0, wirez2, lenLX / 2, lenLY / 2, thickZ1 / 2);
  cathode.SetBoundaryPotential(wireV1);
  geo.AddSolid(&cathode, &Cu);
  

  // Anode plane
  SolidBox anode(0.0, 0.0, -thickZ1/2, lenLX / 2, lenLY / 2, thickZ1 / 2);
  anode.SetLabel("anode");
  anode.SetBoundaryPotential(andplV);
  geo.AddSolid(&anode, &Cu);

  // NeBEM component setup
  
      double tgtElSize = 0.004;  // 0.001
    
  int minEl = 3, maxEl = 5;  // org
  int xcopy = 1, ycopy = 1, zcopy = 0;
  ComponentNeBem3d nebem;
  nebem.SetGeometry(&geo);
  nebem.SetTargetElementSize(tgtElSize);
  nebem.SetMinMaxNumberOfElements(minEl, maxEl);
  nebem.SetPeriodicityX(lenLX);
  nebem.SetPeriodicityY(lenLY);
  nebem.SetPeriodicCopies(xcopy, ycopy, zcopy);
  nebem.UseLUInversion();
  nebem.SetNumberOfThreads(6);
  nebem.Initialise();

 
  ViewGeometry geomView2d;
  geomView2d.SetGeometry(&geo);
  geomView2d.SetArea(xmin, zmin, xmax, zmax);
  geomView2d.SetPlaneXZ();
  TCanvas* c1=new TCanvas("","",600,600);
  geomView2d.SetCanvas(c1);
  
  ViewGeometry geomView2d2;
  geomView2d2.SetGeometry(&geo);
  geomView2d2.SetArea(xmin, zmin, xmax, zmax);
  geomView2d2.SetPlaneXZ();
  TCanvas* c2=new TCanvas("","",600,600);
  geomView2d2.SetCanvas(c2); 

  // Plot device geometry in 3D
  ViewGeometry geomView3d;
  geomView3d.SetGeometry(&geo);
  geomView3d.Plot(); 


  ViewField fieldView;
  fieldView.SetCanvas(c1);
  
  fieldView.SetComponent(&nebem);
  fieldView.SetPlaneXZ();
  // Set the plot limits in the current viewing plane.

  fieldView.SetArea(xmin, zmin, xmax, zmax);
  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, 0.0048, xmax, 0, 0.0048, xf, yf, zf, 60); // 0.0056,
fieldView.PlotFieldLines(xf, yf, zf, true, false); 
geomView2d.Plot2d();
//  c1->Print("NeBem.ps"); 

  // Create the sensor.
  Sensor sensor;
  sensor.AddComponent(&nebem);


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


  ViewDrift driftView;
  constexpr bool plotDrift = true;
  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) { 
//    std::cout << i << "/" << nEvents << "\n";
//    Randomize the initial position. 
      double x0 = 0.0; //-0.5 * wirepitch + RndmUniform() * wirepitch;
      double y0 = 0.0; //-0.5 * wirepitch + RndmUniform() * wirepitch;
//      double z0 = 0.0048; //0.0056; 
      double z0 = 0.0049;
      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.0056, tee = 0.;
      int status = 0;

  int ElectronNumber=0;


//avalMC.DriftElectron(x0, y0, z0, 0.0);
avalMC.AvalancheElectron(x0, y0, z0, 0.0, false);




//avalMC.DisableDiffusion();
//avalMC.DisableAttachment();
//avalMC.ResumeAvalanche(true,false);
//avalMC.GetElectronEndpoint(xee, yee, zee, tee, status);
//avalMC.GetElectronEndpoint(0, x0, y0, z0, t0, xee, yee, zee, tee, status);
const unsigned int np = avalMC.GetNumberOfElectronEndpoints();
    std::cout << np << "/" << nEvents << "\n";
//const unsigned int np2 = avalMC.GetAvalancheSize();   
//    std::cout << np2 << "/" << nEvents << "\n";    
  }
 

  if (plotDrift) {
    TCanvas* c2 = new TCanvas();
 {
      driftView.SetPlane(0, -1, 0, 0, 0, 0);
      driftView.SetArea(xmin, zmin, xmax, zmax); //(-2 * wirepitch, -0.02, 2 * wirepitch, 0.02);
      driftView.SetCanvas(c2);
      constexpr bool twod = true;
      driftView.Plot(twod);
    geomView2d2.Plot2d();
//  c2->Print("NeBem2.ps"); 
    }
  }

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


  
 app.Run(true);
}