#include <cstdlib>
#include <iostream>
#include <fstream>
#include <vector>
#include <math.h>
#include <random>

#include <TCanvas.h>
#include <TROOT.h>
#include <TApplication.h>
#include <TH1D.h>
#include <TGraph.h>
#include <TTree.h>
#include <TFile.h>
#include <TSystem.h>

#include "Garfield/ComponentFieldMap.hh"
#include "Garfield/ComponentGrid.hh"

#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Medium.hh"
#include "Garfield/MediumConductor.hh"

#include "Garfield/Sensor.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/DriftLineRKF.hh"

#include "Garfield/ViewSignal.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ViewDrift.hh"

#include "Garfield/Random.hh"

#include "Garfield/ComponentComsol.hh"

using namespace Garfield;

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

    // This program is to simulate drifts, avalanches and x-ray events in a micromegas chamber

    //Directory of COMSOL field files
    std::string fileloc = "fields_smaller/";

    // Specific field files
    std::string meshfile = "mesh_smaller_new.mphtxt";
    std::string field = "field_smaller_new.txt";
    std::string materials = "materials1.dat";
    std::string gas_file = "build/Ar_90_co2_10_0.0T_293.15K.gas";

    // Directory of output CSV file
    std::string out_file = "Output/MicromegasMain.csv";

    // Specify number of readout strips in single axis
    int numReadoutstrips = 5;


    // Assign directories 
    meshfile = fileloc + meshfile;
    field = fileloc + field;
    materials = fileloc + materials;

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

    // Import field map

    std::cout<<"Beginning to read COMSOL"<<"\n";

    ComponentComsol fm(meshfile, materials, field, "cm");
    fm.EnableConvergenceWarnings(false);

    std::cout<<"Done reading COMSOL"<<"\n";

    // Number of events
    int numEvents = 1;

    // Begin to implement weighting fields
    // Create vector with labels of each weighting fields 
    std::vector<std::string> w_labels;
    // Define x weighting fields
    for (int i = 0; i < numReadoutstrips; ++i) {
        w_labels.push_back("x" + std::to_string(i+1));
        continue;
    }
    // Define y weighting fields
    for (int i = 0; i < numReadoutstrips; ++i) {
        w_labels.push_back("y" + std::to_string(i+1));
        continue;

    }
    std::cout << "Done loading weighting fields, testing weighting field vector labels \n";

    for (int i = 0; i < w_labels.size(); ++i ) {
        std::cout << "Index: " << i << " is weighting label: " << w_labels[i] << "\n"; 
    }

    // Note for future, this can take a while, if you're using weighting fields comment this out
    // load weighting potentials
	for (int i = 0; i < w_labels.size(); ++i) {
		std::cout<<"Reading weighting field "<<w_labels[i]<<"... \n";
        // Indicate name of file containing weighting field
		std::string w_field = fileloc+"field-"+w_labels[i]+".txt";
		
		fm.SetWeightingField(w_field, w_labels[i]);
	}





    // Import gas medium
    MediumMagboltz gas;
    gas.LoadGasFile(gas_file);
    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_CO2+_CO2.txt");
    gas.EnableDrift();

    gas.Initialise(true);

    //assign materials to their respective domains
	const unsigned int nMaterials = fm.GetNumberOfMaterials();
  	for (unsigned int i = 0; i < nMaterials; ++i) {
    	const double eps = fm.GetPermittivity(i);
    	if (eps == 1.) fm.SetMedium(i, &gas);
    }

    // Create sensors for each electrode
    Sensor sensor;
    sensor.AddComponent(&fm);
    // Add electrodes to sensor
    for (int i = 0; i < w_labels.size(); ++i) {
        sensor.AddElectrode(&fm, w_labels[i]);
    }

    // Setup other stuff for signals
    // Want to avalanche ions
    AvalancheMC drift;
    drift.SetSensor(&sensor);
    drift.SetDistanceSteps(2.e-4);

 

    // Calculating Avalanches
    AvalancheMicroscopic aval;
    aval.SetSensor(&sensor);
    aval.UseWeightingPotential(true);

    // Begin plotting drift lines
    // ViewDrift driftView;
    // TCanvas* cDrift = new TCanvas("cDrift", "", 600, 600);
    // driftView.SetCanvas(cDrift);
    // aval.EnablePlotting(&driftView);


    // Begin drift line plotting
    ViewDrift* driftView = nullptr;
	TCanvas* cDrift = nullptr;
	cDrift = new TCanvas("cDrift", "", 600, 600);
	driftView = new ViewDrift();
	driftView->SetCanvas(cDrift);
	aval.EnablePlotting(driftView);
    drift.EnablePlotting(driftView);

    // Begin signal plotting
    ViewSignal* signalView = nullptr;
    TCanvas* cSignal = nullptr;
    cSignal = new TCanvas("cSignal", "", 600, 600);
    signalView = new ViewSignal();
    signalView -> SetCanvas(cSignal);
    signalView -> SetSensor(&sensor);
    sensor.EnableDebugging();


    // Set time binning for signal
    // Note: tmin and tmax provide the times at which you want to record the signal betwee and ntimebins tells you how many bins
    // you want during this interval e.g tmin = 0, tmax = 100, nTimeBins = 100 is the same as recording the signal between 
    // t = 0 and 100 ns with 1 bin per ns
    // For this reason you would probably want to start recording the signal when the electron passes the mesh and actually produces
    // a signal
    const unsigned int nTimeBins = 100;
    const double tmin = 0;
    const double tmax = 100;
    const double tstep = (tmax - tmin) / nTimeBins;
    sensor.SetTimeWindow(tmin, tstep, nTimeBins);

    // Begin plotting signal
    // ViewSignal signalView;
    // TCanvas* cSignal = new TCanvas("cSignal", "", 600, 600);
    // signalView.SetCanvas(cSignal);
    // signalView.SetSensor(&sensor);


    // Assign values types
    int ne, ni;


    // Assign range to randomly generate numbers
    double range_from = 0.15;
    double range_to = 0.25;

    // Create a random number generator for generating large sample sizes
    std::default_random_engine generator;
    std::uniform_real_distribution<double> distribution(range_from, range_to);

    // Produce csv file to write data to
    std::ofstream ava_file(out_file);

    // Build header for csv file
    ava_file << "Testing csv file \n";


    // Begin processing events

    for (int eventInter = 0; eventInter < numEvents; ++eventInter) {
        sensor.ClearSignal(); // reset signal to 0 after each event


        std::cout << "------------ Processing Event: " << eventInter << "/" << numEvents - 1 << " ----------------\n";

        ava_file << eventInter << ",";


        // Assign initial values

        // If random values


        // If specific point
        // double xi = distribution(generator);
        double xi = 0.22;
        // double yi = distribution(generator);
        double yi = 0.22;
        double zi = 0.25;
        double ti = 0;
        double ei = 0;

        ava_file << "xi: " << xi << "," << "yi: " << yi << "," << "zi: " << zi << ",";

        // Send the electron on an avalanche (goodbye electron)

        // std::cout << "xi = " << xi << " yi = " << yi << "\n";
        aval.AvalancheElectron(xi, yi, zi, ti, ei, 0, 0, -1);
        // The above is with the initial electron travelling down (direction given by last 3 coordinates)
        aval.GetAvalancheSize(ne, ni);
        std::cout << "ne = " << ne << " ni = " << ni << "\n";

        // Loop over all electrons in avalanche
        for (const auto& electron : aval.GetElectrons()) {
            // Simulate an ion drift line from the starting point of the electron
            const auto& p0 = electron.path[0];
            drift.DriftIon(p0.x, p0.y, p0.z, p0.t);

        }

        ava_file << "ne: " << ne << "," << "ni: " << ni << "\n";


        // Plot
    //     driftView.SetPlaneXZ();
    //     driftView.SetArea(0.15, 0.14, 0.3, 0.3);
    //     driftView.Plot();
        driftView -> SetPlaneXZ();
	driftView -> SetArea(0, 0, 0, 0.5, 0.5, 0.6);
	driftView -> Plot(true);



        // Plotting signal related tings
        signalView -> PlotSignal("x3");


    }
    ava_file.close();

    app.Run(true);

}