#include <iostream>
#include <fstream>
#include <cmath>
#include "math.h"

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

#include "Garfield/MediumConductor.hh"
#include "Garfield/ViewGeometry.hh"
#include "Garfield/SolidBox.hh"
#include "Garfield/ComponentAnsys123.hh"
#include "Garfield/ComponentComsol.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ViewFEMesh.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/Random.hh"
#include "Garfield/TrackSrim.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/ViewCell.hh"
#include "Garfield/Plotting.hh"
#include <TSystem.h>
#include <cstdlib>
#include <iostream>

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

#include "Garfield/ComponentAnsys123.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ViewFEMesh.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/Random.hh"


using namespace Garfield;
using namespace std;

void ModifyLineData(string fileName, int lineNum, char* lineData)
{
	cout<<"ModifyLineData"<<endl;
	ifstream in;
	in.open(fileName.c_str());
	string strFileData = "";
	int line = 1;
	char tmpLineData[1024] = {0};
	while(in.getline(tmpLineData, sizeof(tmpLineData)))
	{
		if (line == lineNum)
		{
			strFileData += string(lineData);
			strFileData += "\n";
		}
		else
		{
			strFileData += string(tmpLineData);
			strFileData += "\n";
		}
		line++;
	}
	in.close(); 
	//写入文件
	ofstream out;
	out.open(fileName.c_str());
	out.flush();
	out<<strFileData;
	out.close();
}


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

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


	string filename;
	string mypath;
	string commend;

	//
	// ---- 绘图/存储控制 ----
	//
	constexpr bool plotField = true;
	constexpr bool plotDrift = true;
	constexpr bool plotSignal = true;
	constexpr bool savefile = true;

	// Dimension

	double Dim_x = 20;
	double Dim_y = 20;
	
	double Dim_z_down = -25; 
	double Dim_z_up = 25 ;    // cm


	//
	// ---- Track ----
	//
	
	double x0 = 0.5;
	double y0 = 0.5;
	double z0 = 0.5;
	double t0 = 0.;
	double dx = 1.;
	double dy = 1.;
	double dz = 1.;


	double Total_counts;
	double Counts;

  	// 簇信息
	double xc = 0., yc = 0., zc = 0., tc = 0., ec = 0., ekin = 0.;
	int nc = 0, status = 0;

  	// electrons start points and stop points 
	double xe1 = 0., ye1 = 0., ze1 = 0., te1 = 0., ee1 = 0.;
	double xe2 = 0., ye2 = 0., ze2 = 0., te2 = 0., ee2 = 0.;
 	

	string file_pAnode1, file_pAnode2, file_pCathode,file_mesh, file_field;

	//
	// ---- 画板 ----
	//

	TCanvas* ChargeSignal_pAnode1 = nullptr;
	TCanvas* ChargeSignal_pAnode2 = nullptr;
	TCanvas* ChargeSignal_pCathode = nullptr;

	TCanvas* cFieldXY = nullptr;
	TCanvas* cFieldXZ = nullptr;
	TCanvas* cDrift = nullptr;

 	 // Dimension


	Total_counts = 0;
	Counts = 0;


  	// Setup the gas.

	MediumMagboltz gas;
	gas.LoadGasFile("./gasfile/gas.gas");
	const std::string path = std::getenv("GARFIELD_INSTALL");
	gas.LoadIonMobility(path + "/share/Garfield/Data/IonMobility_Ar+_Ar.txt");
	

	//
	// ---- Load the field map ----
	//

	const std::string label_pAnode1 = "readout_pAnode1";
	const std::string label_pAnode2 = "readout_pAnode2";
	const std::string label_pCathode = "readout_pCathode";

	file_pAnode1 =   "./comsolfile/pAnode1.txt";
	file_pAnode2 =   "./comsolfile/pAnode2.txt";
	file_pCathode =   "./comsolfile/pCathode.txt";
	file_field=   "./comsolfile/wfield.txt";
	file_mesh=   "./comsolfile/mesh.mphtxt";

	ModifyLineData(file_pAnode1.c_str(),
	9,
	"% x                       y                        z                        V (V)");

	ModifyLineData(file_pAnode2.c_str(),
	9,
	"% x                       y                        z                        V (V)");

	ModifyLineData(file_pCathode.c_str(),
	9,
	"% x                       y                        z                        V (V)");

	ModifyLineData(file_field.c_str(),
	9,
	"% x                       y                        z                        V (V)");

	//
	// ---- Load the field map ----

	ComponentComsol fm;
	fm.Initialise(file_mesh.c_str(), "./gasfile/dielectrics.dat", file_field.c_str(), "cm");
	fm.SetWeightingField(file_pAnode1.c_str(),label_pAnode1);
	fm.SetWeightingField(file_pAnode2.c_str(),label_pAnode2);
	fm.SetWeightingField(file_pCathode.c_str(),label_pCathode);
	fm.PrintRange();

	// Associate the gas with the corresponding field map material. 
	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);
	}
	fm.PrintMaterials();

	//
	// ---- Signal时间窗 ----
	//	
	constexpr double tMin = 0.;
	constexpr double tMax = 10000.;   // ns
	constexpr double tStep = 1;
	constexpr int nTimeBins = int((tMax-tMin)/tStep);

	// Create the sensor.
	Sensor sensor;
	sensor.Clear();
	
	sensor.AddComponent(&fm);
	sensor.SetArea(-Dim_x, -Dim_y, Dim_z_down,
					Dim_x, Dim_y,  Dim_z_up);

	sensor.AddElectrode(&fm, label_pAnode1);
	sensor.AddElectrode(&fm, label_pAnode2);
	sensor.AddElectrode(&fm, label_pCathode);
	sensor.SetTimeWindow(0., tStep, nTimeBins);


	//
	// ---- 绘制电场 ----
	//

	ViewField* fieldViewXY = new ViewField();
	ViewField* fieldViewXZ = new ViewField();

	if(plotField) 
	{

		cFieldXZ = new TCanvas("cFieldXZ", "", 800, 800);
		fieldViewXZ->SetComponent(&fm);
  		fieldViewXZ->SetPlane(0, -0.1, 0, 0, 0.05, 0);
		fieldViewXZ->SetArea(-Dim_x, Dim_z_down, Dim_x, Dim_z_up);
		fieldViewXZ->SetSensor(&sensor);
		fieldViewXZ->SetNumberOfSamples2d(500,500);
		//fieldViewXZ->SetElectricFieldRange(0,2000);
		fieldViewXZ->SetCanvas(cFieldXZ);
		fieldViewXZ->Plot("e","colz");

		cFieldXY = new TCanvas("cFieldXY", "", 900, 800);
		fieldViewXY->SetComponent(&fm);
  		fieldViewXY->SetPlane(0, 0, -0.1,  0, 0, 0.05);
		fieldViewXY->SetArea(-Dim_x, -Dim_y, Dim_x, Dim_y);
		fieldViewXY->SetSensor(&sensor);
		fieldViewXY->SetNumberOfSamples2d(500,500);
		//fieldViewXY->SetElectricFieldRange(0,2000);
		fieldViewXY->SetCanvas(cFieldXY);
		fieldViewXY->Plot("e","colz");

		//fieldViewXY->PlotWeightingField(label_Anode1,"p","colz");
	}
  

	//
	// ---- Track ----
	//	
	TrackSrim *tr = new TrackSrim();
	tr->SetSensor(&sensor);
	tr->ReadFile("./gasfile/alpha.txt");						// 导入气体文件
	tr->SetKineticEnergy(5.486e6);							    	// 粒子能量5.486MeV
	tr->SetWorkFunction(26580);									// 平均电离能 , 为了计算速度加快，在此放大了1000倍
	tr->SetFanoFactor(0.3);										// 法诺因子

	constexpr double za = 0.9 * (18. / 40.) + 0.1 * (10. / 16.);
	tr->SetAtomicMassNumbers(10. / za, 10.);					// 气体靶的A/Z值
	tr->SetClustersMaximum(100);								// 一个团簇多少电子
	tr->Print();

	//
	// ---- Micro 电子漂移轨迹 ----
	//

	//AvalancheMicroscopic *aval = new AvalancheMicroscopic();
	AvalancheMC*aval = new AvalancheMC();
	//aval->EnableWeightingFieldIntegration();

	aval->SetSensor(&sensor);
	//aval->EnableAttachmentMarkers();
	aval->EnableSignalCalculation(true);
	aval->UseWeightingPotential(true);
	//aval->SetTimeWindow(tMin, tMax);
	//aval->SetElectronTransportCut(0.);
	aval->DisableAvalancheSizeLimit();	
	//aval->EnableDriftLines();

	//
	// ---- 几何、track、电子轨迹 可视化 ----
	//
/*
	//ViewCell cellView;
    TCanvas canvas("cCanvas", "", 600, 600);
    canvas.SetLeftMargin(0.16);
*/
	ViewDrift *driftView = new ViewDrift();
	if (plotDrift) 
	{	
		cDrift = new TCanvas("cDrift", "", 600, 600);

		driftView->SetClusterMarkerSize(0.1);
		driftView->SetCollisionMarkerSize(0.1);
		driftView->SetColourElectrons(kRed);
		driftView->SetColourTracks(kGreen + 3);
		tr->EnablePlotting(driftView);
		aval->EnablePlotting(driftView);
		//driftView->SetCanvas(&canvas);
	}
  
  TRandom rrand;
  double rand;
  const unsigned int nTracks = 30;   // ---  释放的Alpha 粒子个数 --- //
  TLatex label;
  int count = 0;


  for (unsigned int j = 0; j < nTracks; ++j) 
  {
    //sensor.ClearSignal();
	cout << " ------------------------  New Track ------------------------ "<<endl;

	rand = rrand.Uniform(0,360);
	x0 = 0.1*sin(rand*Pi);                  //  Alpha出射位置（随机）
	y0 = 0.1*cos(rand*Pi);

	z0 = -4.99;             // Ａnode2和Cathode之间释放 

	dx = rrand.Uniform(-1,1);     // Alpha正方向出射（随机）
	dy = rrand.Uniform(-1,1);
	dz = rrand.Uniform(0,1);

	tr->NewTrack(x0, y0, z0, t0, dx, dy, dz);
	// loop clusters

	while(tr->GetCluster(xc, yc, zc, tc, nc, ec, ekin))        //获取当前团簇信息
	{
		ee1 = 0.1; 		 //当前cluster的平均电子初始能量

		cout << " ===================================================== "<<endl;
		cout << " Number of Electron: " << nc <<endl;           					//当前cluster的电子个数
		cout << " Position of Clusters: " << xc <<","<< yc <<","<< zc <<endl;		//当前cluster的位置

		
		// loop electrons
		for(int i = 0; i < nc; i++)                //初电子循环
		{
			//int ne = 0, ni = 0;											 
			//aval->GetAvalancheSize(ne, ni);								//雪崩电子数,(不雪崩）
			//driftView->Clear();

				aval->DriftElectron(xc, yc, zc, tc);	//开始计算电子，ee1，初电子各向同性
				const int np = aval->GetNumberOfElectronEndpoints();     	//倍增电子个数，1
				cout << " Multiplicity: " << np << endl;	

				for(int j = 0; j < np; j++)
				{			
					aval->GetElectronEndpoint(j, xe1, ye1, ze1, te1, xe2, ye2, ze2, te2, status);    //倍增电子最终状态
					cout << " Electron start position: " << xe1 <<","<< ye1 <<","<< ze1 <<"\t"<< " Electron stop position: " << xe2 <<","<< ye2 <<","<< ze2 <<endl;
					//eTrack->Fill();
					Total_counts++;
				}
		}
	}
  }
  


	//
	// ---- 几何、track、电子轨迹 可视化 ----
	//

	if (plotDrift) {
	  
    constexpr bool plotMesh = false;
    if (plotMesh) {
      ViewFEMesh* meshView = new ViewFEMesh();
      meshView->SetArea(-Dim_x/20, -Dim_y/20, Dim_z_down,
					     Dim_x/20,  Dim_y/20, Dim_z_up);
				
      meshView->SetCanvas(cDrift);
      meshView->SetComponent(&fm);
      // x-z projection.
      meshView->SetPlane(0, -0.1, 0, 0, 0.05, 0);
      meshView->SetFillMesh(true);

	  driftView->SetColourElectrons(kRed);
	  driftView->SetColourTracks(kGreen + 3);     

      // Set the color.
      meshView->SetViewDrift(driftView);
      meshView->Plot();
    } else {

      driftView->SetPlane(0, -0.1, 0, 0, 0, 0);
      driftView->SetArea(-Dim_x/4., Dim_z_down/4.,
				        Dim_x/4., Dim_z_up/4.);
      driftView->SetCanvas(cDrift);
      constexpr bool twod = true;
      driftView->Plot(twod, twod);
    }
  }


	ofstream signal_pAnode1;
	signal_pAnode1.open(  "SignalpAnode1.txt");

	ofstream signal_pAnode2;
	signal_pAnode2.open(  "SignalpAnode2.txt");

	ofstream signal_pCathode;
	signal_pCathode.open(  "SignalpCathode.txt");

	double Q_Anode1=0;
	double Q_Anode2=0;
	double Q_Cathode=0;

	for(int i = 0; i < nTimeBins; i++)
	{
		Q_Anode1 = (sensor.GetSignal(label_pAnode1, i));
		signal_pAnode1 << i << "\t" << Q_Anode1 << endl;

		Q_Anode2 = (sensor.GetSignal(label_pAnode2, i));
		signal_pAnode2 << i << "\t" << Q_Anode2 << endl;

		Q_Cathode = (sensor.GetSignal(label_pCathode, i));
		signal_pCathode << i << "\t" << Q_Cathode << endl;

	}
	
	//
	// ---- 输出电荷信号可视化 ----
	//

	ViewSignal Signalview_pAnode1;                      //  pAnode1 信号可视化
	ViewSignal Signalview_pAnode2;                      //  pAnode2 信号可视化
	ViewSignal Signalview_pCathode;                    

	if (plotSignal) 
	{
		ChargeSignal_pAnode1 = new TCanvas("Anode1 Signal", "", 600, 600);
		Signalview_pAnode1.SetSensor(&sensor);
		sensor.IntegrateSignal(label_pAnode1);
		Signalview_pAnode1.SetCanvas(ChargeSignal_pAnode1);
		Signalview_pAnode1.SetLabelY("Induced Charge [fC]");
		Signalview_pAnode1.PlotSignal(label_pAnode1,true,false,false,true,false);

		ChargeSignal_pAnode2 = new TCanvas("Anode2 Signal", "", 600, 600);
		Signalview_pAnode2.SetSensor(&sensor);
		sensor.IntegrateSignal(label_pAnode2);
		Signalview_pAnode2.SetCanvas(ChargeSignal_pAnode2);
		Signalview_pAnode2.SetLabelY("Induced Charge [fC]");
		Signalview_pAnode2.PlotSignal(label_pAnode2,true,false,false,true,false);

		ChargeSignal_pCathode = new TCanvas("Cathode Signal", "", 600, 600);
		Signalview_pCathode.SetSensor(&sensor);
		sensor.IntegrateSignal(label_pCathode);
		Signalview_pCathode.SetCanvas(ChargeSignal_pCathode);
		Signalview_pCathode.SetSensor(&sensor);
		Signalview_pCathode.SetLabelY("Induced Charge [fC]");
		Signalview_pCathode.PlotSignal(label_pCathode);

	}


	if(savefile)
	{	
		string rootfile = "/Results.root";
		TFile *hist = new TFile("Results.root", "RECREATE");
		ChargeSignal_pAnode1->Write("Anode1 Signal");
		ChargeSignal_pAnode2->Write("Anode2 Signal");
		ChargeSignal_pCathode->Write("Cathode Signal");
		cFieldXY->Write("FieldXY");
		cFieldXZ->Write("FieldXZ");
		cDrift->Write("Drift");
		hist->Close();
	}

	cout << "=======================" <<endl;
	cout << " The code end! " << endl;

	app.Run(kTRUE);
	return 0;
}