#include <iostream>
#include <fstream>
#include <sstream>

#include <TApplication.h>
#include <TCanvas.h>
#include <TROOT.h>
#include <TFile.h>
#include <TH2F.h>
#include <TH1F.h>

#include "MediumMagboltz.hh"
#include "SolidBox.hh"
#include "GeometrySimple.hh"
#include "ComponentAnalyticField.hh"
#include "Sensor.hh"
#include "TrackHeed.hh"
#include "Plotting.hh"
#include "ViewGeometry.hh"
#include "ViewCell.hh"
#include "Random.hh"
#include "DriftLineRKF.hh"
#include "AvalancheMC.hh"
#include "AvalancheMicroscopic.hh"

using namespace Garfield;
using namespace std;

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

	MediumMagboltz* gas = new MediumMagboltz();
	gas->SetComposition("ar", 80., "co2", 20.);	// Component
	gas->SetTemperature(293.15); // K
	gas->SetPressure(750);	// Torr

	const double width = 2.;	// Gap [cm]
	SolidBox* box = new SolidBox(0., 0., 0., 10., width / 2., 10.);	// y direction	
	GeometrySimple* geo = new GeometrySimple();
	geo->AddSolid(box, gas);
	
	ComponentAnalyticField* comp = new ComponentAnalyticField();
	// add plane
	comp->AddPlaneY(-width / 2 ,0.,"Cathode");
	comp->AddPlaneY(width / 2, 1000. ,"Anode");
	// magnetic field T
	comp->SetMagneticField(0.6, 0., 0.);
	comp->SetGeometry(geo);	
/*	
	// View Cell	
    TCanvas* cellcanvas = new TCanvas();	
	ViewCell* cellView = new ViewCell();// Create a viewer.	
	cellView->SetComponent(comp);// Set the pointer to the component.
    cellView->SetArea(-10., -width/2., -10., 10., width/2., 10);
    cellView->SetCanvas(cellcanvas);	
	cellView->Plot2d();// Make a two-dimensional plot of the cell layout.
    cellcanvas->Update();
  
    ViewDrift* driftview = new ViewDrift();
    driftview->SetArea(-10., -width/2., -10., 10., width/2., 10);
    driftview->SetCanvas(cellcanvas);
*/
	Sensor* sensor = new Sensor();
	sensor->AddComponent(comp);
	
	TrackHeed* track = new TrackHeed();
	track->SetSensor(sensor);
//	track->EnablePlotting(driftview);
	track->SetParticle("p");
	// Set the Kinetic Energy eV.
	const double KineticEnergy = 1.e9;
	track->SetKineticEnergy(KineticEnergy);
	double x0=0., y0=-width/2., z0=0., t0=10.;
	double dx0=0., dy0=1., dz0=0.;
	track->NewTrack(x0, y0, z0, t0, dx0, dy0, dz0);	
/*	
	// particle track
    driftview->Plot(true,false);	//2D
    cellcanvas->Update();	
*/	
	AvalancheMicroscopic* microscopic = new AvalancheMicroscopic();
	microscopic->SetSensor(sensor);
	microscopic->EnableMagneticField();
/*
    microscopic->EnablePlotting(driftview);	// plot drift line
    microscopic->EnableDriftLines();	//Switch on storage of drift lines	
*/	
	TFile *f= new TFile("output.root","RECREATE");
	f->cd();	
	TH1F* eene = new TH1F("eene", "electron end energy", 100, 0, 10);
	TH2F* epos = new TH2F("epos", "electron end position", 100, -5, 5,100,-5,5);	
	
	double xcl, ycl, zcl, tcl, ecl, extra;
	int ncl = 0;
	double ne=0; // primary electron count
	while (track->GetCluster(xcl, ycl, zcl, tcl, ncl, ecl, extra))
	{
		ne += ncl;
		for(int i = 0; i < ncl; i++)
		{
			double x, y, z, t, e, dx, dy, dz;
			track->GetElectron(i,x,y,z,t,e,dx,dy,dz);
			
			// microscopic
			// drift electron
			microscopic->DriftElectron(x,y,z,t,e,dx,dy,dz);
/*			
	        //plot all electron track
            driftview->Plot(true,false);	//2D
            cellcanvas->Update();
*/			
			// get electron end point
			double ex0,ey0,ez0,et0,ee0,ex1,ey1,ez1,et1,ee1;
			int status;
			microscopic->GetElectronEndpoint(0,ex0,ey0,ez0,et0,ee0,ex1,ey1,ez1,et1,ee1,status);
			
			// fill histograms
			eene->Fill(ee1);
			epos->Fill(ex1,ez1);
		}
	}
  
/*
	TCanvas* distr = new TCanvas(800,600);
	distr->Divide(1,2);
	distr->cd(1);
	epos->Draw();
	distr->cd(2);
	eene->Draw();
*/	
	f->Write();
	f->Close();
	cout<<"total electron: "<<ne<<endl;
	
	app.Run(kTRUE);
	
}