#include <iostream>
#include <fstream>
#include <cstdlib>
#include <sys/stat.h>
#include <math.h>

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

#include "Garfield/ComponentNeBem3d.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/MediumConductor.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/SolidBox.hh"
#include "Garfield/SolidWire.hh"
#include "Garfield/ViewGeometry.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewSignal.hh"

using namespace Garfield;

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

    //load the gas medium for the internal space
    MediumMagboltz gas;
    gas.LoadGasFile("ar_85_co2_15_1atm.gas");
    //define the conductor material for the strips and wire
    MediumConductor plane_medium;
    MediumConductor strip_medium;
    MediumConductor wire_medium;

    //Defining the finite plane geometry through a GeometrySimple object
    GeometrySimple geo;

    //geometric parameters
    //----whole cell----
    //half-length of the assembly along z [cm]
    const double hlength = 100.;
    //gap half-length
    const double half_gap = 0.1;    
    //half width of a cell
    const double cell_width = 1.;
    //half thickness of a cell
    const double cell_thick = 0.5*cell_width;

    //----strips----
    //half dimensions
    const double stripWidth = 1*cell_width-half_gap;
    const double stripThick = 35e-4;


    //----wires----
    //signal wire radius [cm]
    const double r_signal_wire=20.e-4/2;
    //field wire radius
    const double r_field_wire=100.e-4/2;

    //----Voltages----
    //strip voltage
    const double vStrip = -1000.;
    //signal wire voltage
    const double v_signal_wire = +1500.;
    //field wire voltage
    const double v_field_wire = -3500.;

    //define the cathode strips
    SolidBox c_strip_up(0.,cell_thick-stripThick,0.,stripWidth,stripThick,hlength);
    SolidBox c_strip_down(0.,-cell_thick+stripThick,0.,stripWidth,stripThick,hlength);
    //define the signal wire
    SolidWire sens_wire(0.,0.,0.,r_signal_wire,hlength);
    //set this wire as the electrode
    sens_wire.SetLabel("s");
    
    //define the field wires: left side
    SolidWire field_wire_u_l(-cell_width,0.,0.,r_field_wire,hlength);
    //define the field wires: right side
    SolidWire field_wire_u_r(+cell_width,0.,0.,r_field_wire,hlength);
    
    //set the boundary potentials for the components
    sens_wire.SetBoundaryPotential(v_signal_wire);
    field_wire_u_l.SetBoundaryPotential(v_field_wire);
    field_wire_u_r.SetBoundaryPotential(v_field_wire);
    c_strip_up.SetBoundaryPotential(vStrip);
    c_strip_down.SetBoundaryPotential(vStrip);

    //add wires strips and planes to the geometry
    geo.AddSolid(&sens_wire,&wire_medium);
    geo.AddSolid(&field_wire_u_l,&wire_medium);
    geo.AddSolid(&field_wire_u_r,&wire_medium);
    geo.AddSolid(&c_strip_up, &strip_medium);
    geo.AddSolid(&c_strip_down, &strip_medium);

    //add the boundary volume filled with gas
    geo.SetMedium(&gas);
    
    //Before moving on one needs to compute the field with neBEM
    ComponentNeBem3d nebem;
    nebem.SetGeometry(&geo);
    nebem.SetTargetElementSize(0.1);
    nebem.SetMinMaxNumberOfElements(20,20);
    nebem.SetNumberOfThreads(8);
    nebem.SetPeriodicityX(2*cell_width);
    nebem.Initialise();
    
    //Use a Sensor object to interface to the transport classes
    Sensor sensor;
    //add the field and electrode components
    sensor.AddComponent(&nebem);
    sensor.AddElectrode(&nebem,"s");

    //setup the binning for the signal calculation (start, bin width[ns],bin number)
    sensor.SetTimeWindow(0.,1.,300);
    sensor.SetArea(-cell_width,-cell_thick,-hlength,cell_width,cell_thick,hlength);

    //Simulating the ionization with Heed
    //first a track object is defined
    std::cout<<"-- Particle track simulation --"<<std::endl;
    TrackHeed track;
    track.SetParticle("muon");
    track.SetEnergy(170e9);
    track.SetSensor(&sensor);
    //compute the track
    track.EnableDebugging();
    
    //compute the electron drift lines
    DriftLineRKF drift;
    drift.SetSensor(&sensor);
    drift.UseWeightingPotential(false);
    drift.EnableAvalanche();
    drift.EnableSignalCalculation();

    //Setting up the viewer of the simulated drift lines
    //create a new canvas
    TCanvas* cD = new TCanvas("cD"," ",1100,500);
    //the ViewDrift object plots the cluster coordinates and drift line points
    ViewDrift driftView;
    //set the object with which to plot and the canvas
    driftView.SetCanvas(cD);
    drift.EnablePlotting(&driftView);
    track.EnablePlotting(&driftView);

    //set the initial position: middle-left side of the cell
    const double x0 = -cell_width/2, y0 = -cell_thick, z0 = 0., t0 = 0.;
    //set the track direction and define the track
    const double dx0 = 0., dy0 = 1., dz0 = 0.;
    track.NewTrack(x0,y0,z0,t0,dx0,dy0,dz0);

    //loop over the clusters
    double xc = 0., yc = 0., zc = 0., tc = 0., ec = 0., extra = 0.;
    int nc =0;
    while(track.GetCluster(xc,yc,zc,tc,nc,ec,extra)){
        //loop over the electrons in the cluster
        for(int k = 0; k<nc; ++k){
            //electron data vars. are defined
            double xe = 0., ye = 0., ze = 0., te = 0., ee = 0.;
            double dx = 0., dy = 0., dz = 0.;
            double enx = 0., eny = 0., enz = 0., ent = 0.;
            int st = 0.;
            track.GetElectron(k, xe, ye, ze, te, ee, dx, dy, dz);
            drift.DriftElectron(xe, ye, ze, te);
            drift.GetEndPoint(enx, eny, enz, ent, st);
            std::cout<<"x: "<<enx<<" y: "<<eny<<" z: "<<enz<<" t: "<<ent<<" stat: "<<st<<std::endl; 
        }
    }

    std::cout<<"Induced charge: "<<sensor.GetTotalInducedCharge("s")<<std::endl;

    //plot the induced current on the wire by the electrons with ViewSignal
    ViewSignal signalView;
    signalView.SetSensor(&sensor);
    signalView.PlotSignal("s");

    //set the object with which to plot and the canvas
    ViewGeometry geoView;
    geoView.SetGeometry(&geo);
    geoView.SetCanvas(cD);
    geoView.Plot2d();
    
    constexpr bool twod = true;
    constexpr bool drawaxis = false;
    driftView.Plot(twod,drawaxis);

    std::cout<<"End of program"<<std::endl;

    app.Run(true);
    return 0;
}