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

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

#include "Garfield/ComponentComsol.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ViewFEMesh.hh"
#include "Garfield/ViewMedium.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/Random.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/DriftLineRKF.hh"


using namespace Garfield;


int main(int argc, char *argv[]){
	TApplication app("app",&argc,argv);
	ComponentComsol fm;
	fm.Initialise("mesh_induced.mphtxt", "dielectric.dat", "Potential.txt", "mm");
	fm.PrintRange();
	
	ViewField fieldView(&fm);
	constexpr bool plotField = true;
	if(plotField){
		fieldView.SetPlane(0, -1, 0, 0, 0, 0);
		fieldView.SetArea(-0.5, -2.5, 0.5, 12.5);
		fieldView.SetVoltageRange(0., -2100.);
		TCanvas* cf = new TCanvas("cf", "", 600, 600);
		cf->SetLeftMargin(0.16);
		fieldView.SetCanvas(cf);
		fieldView.PlotContour();
	}

	MediumMagboltz gas("sf6");
	gas.SetTemperature(293.15);
	gas.SetPressure(760.);
	gas.Initialise(true);
	
	gas.LoadIonMobility("/Users/tianyu/InducedCharge/SeF5-_SeF6.txt");
	fm.SetGas(&gas);
	fm.PrintMaterials();
	

	fm.SetWeightingField("WeightPotential.txt","anode");
	Sensor sensor(&fm);
	sensor.AddElectrode(&fm,"anode");
	const unsigned int nBins = 10000;
	sensor.SetTimeWindow(0, 2000, nBins);
	
	//DriftLineRKF drift;
	AvalancheMC drift;
	drift.SetSensor(&sensor);
	drift.UseWeightingPotential(true);
	drift.SetDistanceSteps(1e-4);
	drift.EnableSignalCalculation();

	
	drift.DriftIon(0,0,9.5,0);
	
	/*
	TrackHeed track(&sensor);
	track.SetParticle("e-");
	track.SetEnergy(2.996e9);

	track.NewTrack();
	*/



	/*
	drift.DriftIon(0,0,0.95,0);

	DriftLineRKF drift1(&sensor);
	ViewMedium mediumView;
	mediumView.SetMedium(&gas);
	mediumView.PlotVelocity("i",'e');
	mediumView.EnableAutoRangeX(0);
	mediumView.SetRangeE(0,1600,1);
	*/	

	ViewDrift driftView;
	constexpr bool plotdrift = true;
	if (plotdrift){
		TCanvas* c2 = new TCanvas("c2", "Drift Path", 600, 600);
		driftView.SetCanvas(c2);
		drift.EnableSignalCalculation();
		drift.EnablePlotting(&driftView);
		driftView.Plot();
	}
	/*
	ViewFEMesh meshView(&fm);
	meshView.SetCanvas(c2);
	meshView.SetColor(0, kRed);
	meshView.SetColor(1, kYellow);
	meshView.SetColor(2, kGray);
	meshView.SetColor(3, kBlue);
	meshView.EnableAxes();
	meshView.SetViewDrift(&driftView);
	meshView.Plot();
	*/


	TCanvas c3("c3","",600,600);
	sensor.IntegrateSignals();
	sensor.PlotSignal("anode", &c3);
	const double q1 = sensor.GetIonSignal("anode", nBins - 1) / ElementaryCharge;
	std::cout<<"Induced charge: " << q1 << " e-\n";	


	TH1F h("", "", 1000, -1111,-1111);
	for(int i = 0; i<1000;i++){
		sensor.ClearSignal();
		double a = gRandom->Uniform(-0.25,0.25);
		double b = gRandom->Uniform(-0.25,0.25);
		drift.DriftIon(a,b,9.5,0);
		sensor.IntegrateSignals();
		h.Fill(sensor.GetIonSignal("anode",nBins-1)/ ElementaryCharge);
	}

	TCanvas c4("c4","his",800,600);
	h.Draw();
	app.Run();
}