#define plot_DCA_calc_C

#include <TGraph.h>
#include <TH1F.h>
#include <TCanvas.h>
#include <TFile.h>

using namespace std;

void plot_DCA_calc(){

	//load the data from a csv into a 2d vector "array" 
	vector<vector<float>> myarr;
	string filename="/icarus/app/users/tboone/icaruscode_dev/test/test_tpcmatch/evtdisp/data.csv";
        ifstream f_1;
        f_1.open(filename.c_str());
        string dummyline_1; getline(f_1,dummyline_1);
        string line_1, field_1; 
        vector<float> myv_1;

        while(getline(f_1,line_1)){
                myv_1.clear();
                stringstream ss_1(line_1);
                while(getline(ss_1,field_1,',')){
                        myv_1.push_back(stof(field_1));
                }//end getting individual columns
                myarr.push_back(myv_1);
        }//end getline(f, line)
        f_1.close();

	//Guide to column values in the array:
	//0=crtTime; 1=start/end?; (2,3,4)=CRT Hit position (x,y,z); (5,6,7)=TPC track start(x,y,z); (8,9,10)=TPC Track end(x,y,z)
	//Transfer CRT hits to arrays so that they can be plotted on the TGraph
	float crtTime[myarr.size()], startend[myarr.size()]; 
	float crthit_x[myarr.size()], crthit_y[myarr.size()], crthit_z[myarr.size()];
	float tpcstart_x[myarr.size()], tpcstart_y[myarr.size()], tpcstart_z[myarr.size()];
	float tpcend_x[myarr.size()], tpcend_y[myarr.size()], tpcend_z[myarr.size()];

//	TLine tpctracks_xy[myarr.size()], tpctracks_xz[myarr.size()], tpctracks_yz[myarr.size()];	

	cout << "\nFill arrays:\n";

	for(int i=0; i<myarr.size(); i++){

		crtTime[i]=myarr[i][0];
		startend[i]=myarr[i][1];

		crthit_x[i]=myarr[i][2];
		crthit_y[i]=myarr[i][3];
		crthit_z[i]=myarr[i][4];

		tpcstart_x[i]=myarr[i][5];
		tpcstart_y[i]=myarr[i][6];
		tpcstart_z[i]=myarr[i][7];
		
		tpcend_x[i]=myarr[i][8];
		tpcend_y[i]=myarr[i][9];
		tpcend_z[i]=myarr[i][10];

//		tpctracks_xy[i]= new TLine(tpcstart_x[i], tpcstart_y[i], tpcend_x[i], tpcend_y[i]);
//		tpctracks_xz[i]= new TLine(tpcstart_x[i], tpcstart_z[i], tpcend_x[i], tpcend_z[i]);
//		tpctracks_yz[i]= new TLine(tpcstart_y[i], tpcstart_z[i], tpcend_y[i], tpcend_z[i]);

	}//end loop filling TVectors

	//create TGraphs to hold the three 2D planes we are going to project onto, just draw CRT hits for now
	TGraph *gr_xy = new TGraph(myarr.size(), crthit_x, crthit_y); gr_xy->SetMarkerStyle(8); gr_xy->SetMarkerColor(4);// gr_xy->SetDrawOption("AP");
	TGraph *gr_xz = new TGraph(myarr.size(), crthit_x, crthit_z); gr_xz->SetMarkerStyle(8); gr_xz->SetMarkerColor(4);// gr_xz->SetDrawOption("AP");
	TGraph *gr_yz = new TGraph(myarr.size(), crthit_y, crthit_z); gr_yz->SetMarkerStyle(8); gr_yz->SetMarkerColor(4);// gr_yz->SetDrawOption("AP");

	gr_xy->GetXaxis()->SetRangeUser(-600,600); gr_xy->GetYaxis()->SetRangeUser(-400,400);
	gr_xz->GetXaxis()->SetRangeUser(-600,600); gr_xz->GetYaxis()->SetRangeUser(-1300,1300);
	gr_yz->GetXaxis()->SetRangeUser(-400,400); gr_yz->GetYaxis()->SetRangeUser(-1300,1300);

	gr_xy->SetTitle("XY Plane");
	gr_xy->GetHistogram()->GetXaxis()->SetTitle("x");
	gr_xy->GetHistogram()->GetYaxis()->SetTitle("y");

	gr_xz->SetTitle("XZ Plane");
	gr_xz->GetHistogram()->GetXaxis()->SetTitle("x");
	gr_xz->GetHistogram()->GetYaxis()->SetTitle("z");

	gr_yz->SetTitle("YZ Plane");
	gr_yz->GetHistogram()->GetXaxis()->SetTitle("y");
	gr_yz->GetHistogram()->GetYaxis()->SetTitle("z");

	TCanvas *c1 = new TCanvas(); c1->SetGrid();
	gr_xy->Draw("AP");
	c1->SaveAs("xy_test.png");

	TCanvas *c2 = new TCanvas(); c2->SetGrid();
	gr_xz->Draw("AP");
	c2->SaveAs("xz_test.png");
	
	TCanvas *c3 = new TCanvas(); c3->SetGrid();
	gr_yz->Draw("AP");
	c3->SaveAs("yz_test.png");

	//Create plots with TPC tracks on top
	for(int i=0; i<myarr.size(); i++){

		//debug line for now
		if(i>0) continue;

		//Build mini-arrays to create lines for TPC tracks
		float tpcsmall_x[2], tpcsmall_y[2], tpcsmall_z[2];
		tpcsmall_x[0] = tpcstart_x[i]; tpcsmall_x[1] = tpcend_x[i];
		tpcsmall_y[0] = tpcstart_y[i]; tpcsmall_y[1] = tpcend_y[i];
		tpcsmall_z[0] = tpcstart_z[i]; tpcsmall_z[1] = tpcend_z[i];

		//Plots containing the points from the track reconstruction
		TGraph *gr_tpcsmall_xy = new TGraph(2,tpcsmall_x,tpcsmall_y); gr_tpcsmall_xy->SetMarkerStyle(21); gr_tpcsmall_xy->SetMarkerColor(2);
		TGraph *gr_tpcsmall_xz = new TGraph(2,tpcsmall_x,tpcsmall_z); gr_tpcsmall_xz->SetMarkerStyle(21); gr_tpcsmall_xz->SetMarkerColor(2);
		TGraph *gr_tpcsmall_yz = new TGraph(2,tpcsmall_y,tpcsmall_z); gr_tpcsmall_yz->SetMarkerStyle(21); gr_tpcsmall_yz->SetMarkerColor(2);
		
		//Extend line so we can actually see it
		float slope_xy, slope_xz, slope_yz, intercept_xy, intercept_xz, intercept_yz;

		slope_xy = (tpcend_y[i]-tpcstart_y[i])/(tpcend_x[i]-tpcstart_x[i]);
		intercept_xy = tpcstart_y[i]-(slope_xy * tpcstart_x[i]);

		slope_xz = (tpcend_z[i]-tpcstart_z[i])/(tpcend_x[i]-tpcstart_x[i]);
		intercept_xz = tpcstart_z[i]-(slope_xz * tpcstart_x[i]);

		slope_yz = (tpcend_z[i]-tpcstart_z[i])/(tpcend_y[i]-tpcstart_y[i]);
		intercept_yz = tpcstart_z[i]-(slope_yz * tpcstart_y[i]);

		//Extend the lines so that they're not short little nubs on the graph
		int numpts = 261;
		float tpc_x[numpts], tpc_y[numpts], tpc_z[numpts];
		for(int j=0; j<numpts; j++){

			//set initial point, we will travel along z and fill in points using formulae from the geometrical info
			if(j==0) tpc_z[j] = -1300;
			tpc_x[j] = (tpc_z[j] - intercept_xz)/slope_xz;
			tpc_y[j] = (tpc_z[j] - intercept_yz)/slope_yz;
			if(j!=numpts-1) tpc_z[j+1] = tpc_z[j]+10;

		}//end loop making larger lines

		//Make plots with the extended lines
		TGraph *gr_tpcline_xy = new TGraph(numpts,tpc_x,tpc_y);/* gr_tpcline_xy->SetMarkerStyle(1);*/ gr_tpcline_xy->SetLineColor(1); gr_tpcline_xy->SetLineWidth(2);
		TGraph *gr_tpcline_xz = new TGraph(numpts,tpc_x,tpc_z);/* gr_tpcline_xz->SetMarkerStyle(1);*/ gr_tpcline_xz->SetLineColor(1); gr_tpcline_xz->SetLineWidth(2);
		TGraph *gr_tpcline_yz = new TGraph(numpts,tpc_y,tpc_z);/* gr_tpcline_yz->SetMarkerStyle(1);*/ gr_tpcline_yz->SetLineColor(1); gr_tpcline_yz->SetLineWidth(2);


		gr_tpcline_xy->GetXaxis()->SetLimits(-600.,600.); gr_tpcline_xy->GetYaxis()->SetLimits(-400.,400.);
		gr_tpcline_xz->GetXaxis()->SetLimits(-600.,600.); gr_tpcline_xz->GetYaxis()->SetLimits(-1300.,1300.);
		gr_tpcline_yz->GetXaxis()->SetLimits(-400.,400.); gr_tpcline_yz->GetYaxis()->SetLimits(-1300.,1300.);
		gr_tpcline_xy->GetXaxis()->SetRangeUser(-600.,600.); gr_tpcline_xy->GetYaxis()->SetRangeUser(-400.,400.);
		gr_tpcline_xz->GetXaxis()->SetRangeUser(-600.,600.); gr_tpcline_xz->GetYaxis()->SetRangeUser(-1300.,1300.);
		gr_tpcline_yz->GetXaxis()->SetRangeUser(-400.,400.); gr_tpcline_yz->GetYaxis()->SetRangeUser(-1300.,1300.);


		TCanvas *cxy = new TCanvas();
		TCanvas *cxz = new TCanvas();
		TCanvas *cyz = new TCanvas();

		cxy->cd(); cxy->SetGrid();
		TMultiGraph *mgr_xy = new TMultiGraph(); mgr_xy->SetTitle("XY with TPC track");
		mgr_xy->GetXaxis()->SetTitle("X (cm)");
		mgr_xy->GetYaxis()->SetTitle("Y (cm)");
		mgr_xy->Add(gr_xy,"AP");
		mgr_xy->Add(gr_tpcsmall_xy,"AP");
		mgr_xy->Add(gr_tpcline_xy,"AL");
		mgr_xy->GetXaxis()->SetRangeUser(-600.,600.);	mgr_xy->GetXaxis()->SetLimits(-600.,600.);
		mgr_xy->GetYaxis()->SetRangeUser(-400.,400.);	mgr_xy->GetYaxis()->SetLimits(-400.,400.);
		mgr_xy->Draw("ALP");
		cxy->Update();

		cxz->cd(); cxz->SetGrid();
		TMultiGraph *mgr_xz = new TMultiGraph(); mgr_xz->SetTitle("XZ with TPC track");
		mgr_xz->GetXaxis()->SetTitle("X (cm)");
		mgr_xz->GetYaxis()->SetTitle("Z (cm)");
		mgr_xz->Add(gr_xz,"AP");
		mgr_xz->Add(gr_tpcsmall_xz,"AP");
		mgr_xz->Add(gr_tpcline_xz,"AL");
		mgr_xz->GetXaxis()->SetRangeUser(-600.,600.);	mgr_xz->GetXaxis()->SetLimits(-600.,600.);
		mgr_xz->GetYaxis()->SetRangeUser(-1300.,1300.);	mgr_xz->GetXaxis()->SetLimits(-1300.,1300.);
		mgr_xz->Draw("ALP");
		cxz->Update();

		cyz->cd(); cyz->SetGrid();
		TMultiGraph *mgr_yz = new TMultiGraph(); mgr_yz->SetTitle("YZ with TPC track");
		mgr_yz->GetXaxis()->SetTitle("Y (cm)");
		mgr_yz->GetYaxis()->SetTitle("Z (cm)");
		mgr_yz->Add(gr_yz,"AP");
		mgr_yz->Add(gr_tpcsmall_yz,"AP");
		mgr_yz->Add(gr_tpcline_yz,"AL");
		mgr_yz->GetXaxis()->SetRangeUser(-400.,400.);	mgr_yz->GetXaxis()->SetLimits(-400.,400.);
		mgr_yz->GetYaxis()->SetRangeUser(-1300.,1300.);	mgr_yz->GetXaxis()->SetLimits(-1300.,1300.);
		mgr_yz->Draw("ALP");
		cyz->Update();

		
	}//end loop plotting TPC tracks

}