// This file was derived from the snake_macro.C file originally
// written by Micah Veilluex. That macro read data from the original
// ASCII scan file and created a couple of TTree's on the fly. 
//
// This version produces the plots, but uses the TTree's created
// by the snake_mktree program.
//
// May 7, 2009  David lawrence

//----------------
// fdc_3Dplots
//----------------
int fdc_3Dplots(const char *fname="plane004 000deg 2008_11_06 09_26_26.root")
{
	gROOT->Reset();
	gStyle->SetPalette(1,NULL);
	gStyle->SetPadRightMargin(0.125);
  
	// Setting some style properties of plot, drawing it, and adding title/labels
	gStyle->SetMarkerStyle(8);
	gStyle->SetMarkerSize(.75);
	gStyle->SetMarkerColor(2);
	gStyle->SetPalette(1);
	gStyle->SetTitleXSize(.03);
	gStyle->SetTitleYSize(.03);
	gStyle->SetTitleXOffset(1.6);
	gStyle->SetTitleYOffset(1.6);
	gStyle->AxisChoice("X");
	gStyle->SetLabelSize(.03);
	gStyle->AxisChoice("Y");
	gStyle->SetLabelSize(.03);
	gStyle->AxisChoice("Z");
	gStyle->SetLabelSize(.03);
	gStyle->SetMarkerSize(.8);
	gStyle->SetMarkerColor(2);
	gStyle->SetMarkerStyle(8);
	gStyle->SetOptStat(0);

	// Open ROOT file and get pointer to scan TTree
	TFile *f = new TFile(fname);
	TTree *scan = (TTree*)gROOT->FindObject("scan");
	
	// Create a canvas to draw in later
	TCanvas *c1 = new TCanvas("c1", "", 1100, 850);
	c1->SetTicks();


	// Define output directory and make sure it exists
	char imagesdir[512];
	sprintf(imagesdir,"images/%s",fname);
	imagesdir[strlen(imagesdir)-5] = 0;
	char cmd[256];
	sprintf(cmd, "mkdir -p \"%s\"", imagesdir);
	system(cmd);

	// Some constants used to calculate bin sizes etc.
	const int keep = 844;
	const int numPasses = 42;  // number of vertical passes done during entire scan
	const double y_dist = 1.18;  // total scan distance (meters)
	double y_datapoints_approx = 867;  // number of datapoints per vertical scan
	const double datapoint_sep = y_dist/y_datapoints_approx;  // total vert distance scanned (meters) / datapoints per vert scan 
	const double scan_sep = .0257;// scan separation distance in meters (0.0254 m = 1 inch)
	const double x_dist = numPasses*scan_sep;  // total horizontal distance between scans
	const double ybin_dimension = .02;  // y-dimension of bins is 20mm
	const int xbins=numPasses-1;  // number of xbins
	const int ybins= (int)(y_dist/ybin_dimension)-1;  // number of bins in each dimension

	// flatness_plot_3D
	char title[512];
	char filename[512];

	sprintf(title, "Flatness plot for: %s", fname);
	TProfile2D *flatness_plot_3D = new TProfile2D("flatness_plot_3D", title, xbins, 0, x_dist, ybins, 0, y_dist);
	flatness_plot_3D->SetYTitle("X Position (meters)");
	flatness_plot_3D->SetXTitle("Y Position (meters)");
	flatness_plot_3D->SetZTitle("Depth (mm)");
	flatness_plot_3D->GetXaxis()->SetLabelSize(.03);
	flatness_plot_3D->GetYaxis()->SetLabelSize(.03);
	flatness_plot_3D->GetZaxis()->SetTitleSize(.03);
	flatness_plot_3D->GetZaxis()->SetLabelSize(.03);
	flatness_plot_3D->GetZaxis()->SetTitleOffset(1.6);
	c1->SetPhi(-60);
	scan->Draw("depth:y/100.0:x/100.0>>flatness_plot_3D", "depth != 0.0 && depth != 6.35", "prof, surf1");
	sprintf(filename, "%s/flatness_plot_3d.gif", imagesdir);
	c1->SaveAs(filename);
	sprintf(filename, "%s/flatness_plot_3d.pdf", imagesdir);
	c1->SaveAs(filename);

	///////////// code to plot contour with datapoints rather than surface //////////////
	scan->Draw("depth:y/100.0:x/100.0>>htemp", "depth != 0.0 && depth != 6.35");
	htemp->SetYTitle("X Position (meters)");
	htemp->SetXTitle("Y Position (meters)");
	htemp->SetZTitle("Depth (mm)");
	htemp->SetMarkerStyle(6);
	htemp->SetMarkerSize(0.5);
	htemp->SetMarkerColor(kRed);
	htemp->Draw();
	sprintf(filename, "%s/flatness_plot_3d_LINES.gif", imagesdir);
	c1->SaveAs(filename);
	sprintf(filename, "%s/flatness_plot_3d_LINES.pdf", imagesdir);
	c1->SaveAs(filename);

	// Draw RMS as a 3D plot.
	// Explicitly set bins with an RMS greater than 50 to zero so
	// they don't obscure the view
	c1->SetTheta(65);
	c1->SetPhi(-60);
	TH2D *depth_rms = (TH2D*)gROOT->FindObject("depth_rms");
	TH2D *depth_rms_trimmed = (TH2D*)depth_rms->Clone("depth_rms_trimmed");
	for(int ibin=1; ibin<=depth_rms_trimmed->GetNbinsX(); ibin++){
		for(int jbin=1; jbin<=depth_rms_trimmed->GetNbinsY(); jbin++){
			if(depth_rms_trimmed->GetBinContent(ibin, jbin)>25.0)depth_rms_trimmed->SetBinContent(ibin, jbin, 0.0);
		}
	}
	depth_rms_trimmed->GetZaxis()->SetRangeUser(0.0, 30.0);
	depth_rms_trimmed->Draw("surf4");
	sprintf(filename, "%s/RMS_plot_3d.gif", imagesdir);
	c1->SaveAs(filename);
	sprintf(filename, "%s/RMS_plot_3d.pdf", imagesdir);
	c1->SaveAs(filename);

#if 0


  TCanvas *c3 = new TCanvas("c3","x-pos vs. RMS", 1100, 850);
  c3->cd();
  c3->SetTicks(1,1);
  tree_err->Draw("error.RMS:error.err_xpos");
  sprintf(title, "RMS of flatness plot bins for cathode plane %s, oriented %s degrees clockwise of vertical", &cathode, &orient);
  htemp->SetTitle(title);
  htemp->SetXTitle("Y Position (meters)");
  htemp->SetYTitle("RMS of Depth (microns)");
  htemp->GetXaxis()->SetLabelSize(.03);
  htemp->GetYaxis()->SetLabelSize(.03);
  sprintf(filename, "%s/RMS_plot_XvRMS.jpg", imagesdirname);
  c3->SaveAs(filename);
  sprintf(filename, "%s/RMS_plot_XvRMS.root", imagesdirname);
  c3->SaveAs(filename);

  TCanvas *c4 = new TCanvas("c4","y-pos vs. RMS", 1100, 850);
  c4->cd();
  c4->SetTicks(1,1);
  tree_err->Draw("error.RMS:error.err_ypos");
  sprintf(title, "RMS of flatness plot bins for cathode plane %s, oriented %s degrees clockwise of vertical", &cathode, &orient);
  htemp->SetTitle(title);
  htemp->SetXTitle("X Position (meters)");
  htemp->SetYTitle("RMS of Depth (microns)");
  htemp->GetXaxis()->SetLabelSize(.03);
  htemp->GetYaxis()->SetLabelSize(.03);
  sprintf(filename, "%s/RMS_plot_YvRMS.jpg", imagesdirname);
  c4->SaveAs(filename);
  sprintf(filename, "%s/RMS_plot_YvRMS.root", imagesdirname);
  c4->SaveAs(filename);

  TCanvas *c5 = new TCanvas("c5","x-pos vs. y-pos vs. max delta", 1100,850);
  c5->cd();
  c5->SetTicks(1,1);
  c5->SetPhi(-60);  
  tree_err->Draw("error.binheight:error.err_ypos:error.err_xpos");
  sprintf(title, "Max #Delta for cathode plane %s, oriented %s degrees clockwise of vertical", &cathode, &orient);
  htemp->SetTitle(title);
  htemp->SetXTitle("Y Position (meters)");
  htemp->SetYTitle("X Position (meters)");
  htemp->SetZTitle("Max #Delta (microns)");
  htemp->GetXaxis()->SetLabelSize(.03);
  htemp->GetYaxis()->SetLabelSize(.03);
  htemp->GetZaxis()->SetTitleSize(.03);
  htemp->GetZaxis()->SetLabelSize(.03);
  htemp->GetZaxis()->SetTitleOffset(1.6);
  sprintf(filename, "%s/maxDelta_plot_3d.jpg", imagesdirname);
  c5->SaveAs(filename);
  sprintf(filename, "%s/maxDelta_plot_3d.root", imagesdirname);
  c5->SaveAs(filename);

  TCanvas *c6 = new TCanvas("c6","x-pos vs. max delta", 1100,850);
  c6->cd();
  c6->SetTicks(1,1);
  tree_err->Draw("error.binheight:error.err_xpos");
  sprintf(title, "Max #Delta for cathode plane %s, oriented %s degrees clockwise of vertical", &cathode, &orient);
  htemp->SetTitle(title);
  htemp->SetXTitle("Y Position (meters)");
  htemp->SetYTitle("Max #Delta (microns)");
  htemp->GetXaxis()->SetLabelSize(.03);
  htemp->GetYaxis()->SetLabelSize(.03);
  sprintf(filename, "%s/maxDelta_plot_XvHeight.jpg", imagesdirname);
  c6->SaveAs(filename);
  sprintf(filename, "%s/maxDelta_plot_XvHeight.root", imagesdirname);
  c6->SaveAs(filename);

  TCanvas *c7 = new TCanvas("c7","y-pos vs. max delta", 1100,850);
  c7->cd();
  c7->SetTicks(1,1);
  tree_err->Draw("error.binheight:error.err_ypos");
  sprintf(title, "Max #Delta for cathode plane %s, oriented %s degrees clockwise of vertical", &cathode, &orient);
  htemp->SetTitle(title);
  htemp->SetXTitle("X Position (meters)");
  htemp->SetYTitle("Max #Delta (microns)");
  htemp->GetXaxis()->SetLabelSize(.03);
  htemp->GetYaxis()->SetLabelSize(.03);
  sprintf(filename, "%s/maxDelta_plot_YvHeight.jpg", imagesdirname);
  c7->SaveAs(filename);
  sprintf(filename, "%s/maxDelta_plot_YvHeight.root", imagesdirname);
  c7->SaveAs(filename);

  hfile->Write();
  //delete hfile;
#endif
  return 0;
}