// snake_macro.cc is a root macro that generates several plots describing the shape of a
// surface.  An input file must be available that lists depths (doubles), broken by lines
// of the form "vertical pass 'n' complete" (indicating the end of a line of continuous
// data.  This data should follow a snaking pattern where every other line of data is in
// backwards order.
//
// The first plot is a simple contour of the surface (x-position vs. y-position vs. depth).
// Then two plots are generated that describe the amount of local topographical variation
// of the surface.  Each line of data is broken up into several bins, and the root mean
// square of each bin is plotted (x-pos vs. y-pos vs. RMS).  The greatest difference in
// height of a bin is also plotted (x-pos vs. y-pos vs. (max_bin_depth - min_bin_depth). 
//
// This code was written to plot data that describes the surface of a cathode plane, so
// in order to have nicely named plots, the filename of the datafile must be in the
// following format (date time plane# orientation):
// plane### ###deg yyyy-mm-dd hh-mm-ss.txt

#include <iostream>
#include <sstream>
#include <string>
#include <TMath.h>
#include <TRandom.h>
#include <TStyle.h>
#include <TCanvas.h>

// this is used to convert strings from input file to doubles:
double strtodouble(const string& what)
{
  istringstream instr(what);
  double val;
  instr >> val;
  return val;
}

int snake_macro(char info[]) {
  gROOT->Reset();
  gStyle->SetPalette(1,NULL);
 
  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

  // These variables may need to be adjusted depending on scan and when new system implemented
  // Get data from file
//  char info[] = "plane999 999deg 2008_10_26 20_28_06.txt";  // will be passed in by snakescanner, contains file and scan info
  //char info[] = "0000_00_00 00_0_00 plane0 00deg";
  char datafile[512];
  string info_s = info;
  sprintf(datafile, "/home/laser/FLATNESS/flatness_data/%s", info);  // name of datafile

  // Getting cathode plane info from datafile's filename
  size_t found = info_s.find("plane");
  char cathode[128], orient[128];
  sprintf(cathode, "%c%c%c", info[(int)found+5], info[(int)found+6], info[(int)found+7]);
  found = info_s.find("deg");
  sprintf(orient, "%c%c%c", info[(int)found-3], info[(int)found-2], info[(int)found-1]);
  cout << "<cathode plane " << cathode << "> <" << orient << " degrees clockwise of vertical>" << endl;

  // Creating data structure for tree
  double xpos[1000000], ypos[1000000], depth[1000000];

  // opening file and creating tree
  FILE *flat_data = fopen(datafile, "r");
  TFile *flat_data_f = new TFile("flat_data.root", "RECREATE");

  // declaring some variables
  char line[512];
  double temp[5000];
  size_t endpass; // determines when done adding data from a particular line of data
  double *bin_data = new double[ybins];
  int *y_datapoints = new int[numPasses];
  int *y_row_datapoints = new int[numPasses];
  // entering data from input file into arrays
  cout << "Getting data from file: " << datafile << endl;
  int x=0, y=0, z=0;
  while ( fgets(&line, 512, flat_data) ) {
    string line_s = line;
    endpass = line_s.find("horizontal");
    // checking to see if done with current line of data
    if (endpass != string::npos) {
      // scans go both up and down, so for half of them we must reverse order of data
      if (x%2 == 0) { 
	for (int j=0; j<y; j++) {
	  xpos[z-y+j] = (double)(x*scan_sep);
	  ypos[z-y+j] = (double)(j*datapoint_sep);
	  depth[z-y+j] = temp[y-j];  // selecting depth data in reverse order
	}
      }
      cout << "Filled tree with data from pass " << x+1 << ".  " << y << " data points total." << endl;
      y_datapoints[x] = z;
      y_row_datapoints[x] = y;
      x++;
      y=0;
    }
    else {
      if (x%2 == 1 && y < y_datapoints_approx) {  // data from these scans gets added to tree as is
	xpos[z]=(double)(x*scan_sep);
	ypos[z]=(double)(y*datapoint_sep);
	depth[z] = strtodouble(line)*1000; // converting to microns here
      }
      else if (x%2 == 0 && y < y_datapoints_approx) temp[y] = strtodouble(line)*1000; // order of this data will be reversed later (converted to microns)
      y++;
      z++;
    }
  }

  // creating data structures for trees
  double *err_xpos, *err_ypos, *cont_xpos, *cont_ypos, *cont_depth, *RMS, *binheight;
  struct error_t {
    Double_t err_xpos;
    Double_t err_ypos;
    Double_t RMS;
    Double_t binheight;
  };
  error_t error;

  struct cont_t {
    Double_t cont_xpos;
    Double_t cont_ypos;
    Double_t depth;
  };
  cont_t cont;

  // declaring trees
  TTree *tree_cont = new TTree("tree_cont", "flatness data from ascii file");
  tree_cont->Branch("cont", &cont, "cont.cont_xpos/D:cont.cont_ypos/D:cont.depth/D");
  TTree *tree_err = new TTree("tree_err", "Bin error data from flatness plot");
  tree_err->Branch("error", &error, "error.err_xpos/D:error.err_ypos/D:error.RMS/D:error.binheight/D");
  
  //dividing data into bins and calculating RMS and bin height differences for each
  int bin_index=0, last_bin_end=ybins, vert_pass=0, iter=0, keep_data;
  double sum=0, sum2=0, bin_ave, min, max;
  for (int i=0; i<z; i++) {
    cont.cont_xpos = xpos[i];
    if (vert_pass%2==0) cont.cont_ypos = ypos[i] - (y_row_datapoints[vert_pass]-keep)*datapoint_sep;
    else cont.cont_ypos = ypos[i] +(y_row_datapoints[vert_pass]-keep)*datapoint_sep;
    cont.depth = depth[i];
    tree_cont->Fill();
    if (bin_index <= ybins) {
      bin_ave += depth[i];
      bin_data[bin_index] = depth[i];
      bin_index++;
    }
    if (bin_index == ybins || i == y_datapoints[vert_pass] ) { // end bin when we reach end of bin size or end of pass
      if (i == y_datapoints[vert_pass]) {
        vert_pass++;
      }
      bin_ave = bin_ave/bin_index;
      min = bin_data[0];
      max = bin_data[0];
      for (int j=0; j<bin_index; j++) {
	sum2 += (bin_data[j]-bin_ave) * (bin_data[j]-bin_ave);
	if (bin_data[j] < min) { min = bin_data[j]; }
	if (bin_data[j] > max) { max = bin_data[j]; }
      }
      error.binheight = max - min;
      error.RMS = sqrt(sum2/bin_index);
      error.err_xpos = xpos[i];
      if (vert_pass%2==0) error.err_ypos = ypos[i] - (y_row_datapoints[vert_pass]-keep)*datapoint_sep;
      else error.err_ypos = ypos[i] + (y_row_datapoints[vert_pass]-keep)*datapoint_sep;
      tree_err->Fill();
      bin_index = 0;
      sum =0;
      sum2 = 0;
    }
      iter++;
  }
  
  // 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);
  /////////////////////////
  // creating plots
  /////////////////////////
  char filename[512];
//  sprintf(filename, "%s %s cathode_%s_oriented_%s_deg_flatness_plot_3d.root", imagesdirname);
  char make_path[512];
  char imagesdirname[512];
  sprintf(imagesdirname, "images/%s", info);
  imagesdirname[strlen(imagesdirname)-4]=0;
  sprintf(make_path, "mkdir -p \"%s\"", imagesdirname);
  system(make_path);
  TFile *hfile = gROOT->FindObject(filename); if (hfile) hfile->Close();
  hfile = new TFile(filename,"RECREATE","Demo ROOT file with histograms");

  TCanvas *c1 = new TCanvas("c1","x-pos vs. y-pos vs. depth", 1100, 850);
  c1->SetTicks(1,1);
  tree_cont->UseCurrentStyle();
  TProfile2D *hist = new TProfile2D("hist", "Flatness Plot", xbins, 0, x_dist, ybins, 0, y_dist);
  char title[512];
  sprintf(title, "Flatness plot for cathode plane %s, oriented %s degrees clockwise of vertical", &cathode, &orient);
  hist->SetTitle(title);
  hist->SetYTitle("X Position (meters)");
  hist->SetXTitle("Y Position (meters)");
  hist->SetZTitle("Depth (microns)");
  hist->GetXaxis()->SetLabelSize(.03);
  hist->GetYaxis()->SetLabelSize(.03);
  hist->GetZaxis()->SetTitleSize(.03);
  hist->GetZaxis()->SetLabelSize(.03);
  hist->GetZaxis()->SetTitleOffset(1.6);
  c1->SetPhi(-60);
  tree_cont->Draw("cont.depth:cont.cont_ypos:cont.cont_xpos>>hist", "cont.depth != 0.0 && cont.depth != 6.35", "prof, surf1");
  sprintf(filename, "%s/flatness_plot_3d.jpg", imagesdirname);
  c1->SaveAs(filename);
  sprintf(filename, "%s/flatness_plot_3d.root", imagesdirname);
  c1->SaveAs(filename);

  ///////////// code to plot contour with datapoints rather than surface ////////////// 
  TCanvas *c1 = new TCanvas("c1","x-pos vs. y-pos vs. depth", 1100, 850);
  c1->SetTicks(1,1);
  tree_cont->UseCurrentStyle();
  char title[512];
  sprintf(title, "Flatness plot for cathode plane %s, oriented %s degrees clockwise of vertical", &cathode, &orient);
  c1->SetPhi(-60);
  tree_cont->Draw("cont.depth:cont.cont_ypos:cont.cont_xpos>>htemp", "cont.depth != 0.0 && cont.depth != 6.35");
  //TH3D *htemp = (TH3D*)gROOT->FindObject("htemp");
  htemp->SetTitle(title);
  htemp->SetYTitle("X Position (meters)");
  htemp->SetXTitle("Y Position (meters)");
  htemp->SetZTitle("Depth (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/flatness_plot_3d_LINES.jpg", imagesdirname);
  c1->SaveAs(filename);
  sprintf(filename, "%s/flatness_plot_3d_LINES.root", imagesdirname);
  c1->SaveAs(filename);
  

  TCanvas *c2 = new TCanvas("c2","x-pos vs. y-pos vs. RMS", 1100, 850);
  c2->cd();
  c2->SetTicks(1,1);
  tree_err->UseCurrentStyle();
  c2->SetPhi(-60);
  tree_err->Draw("error.RMS:error.err_ypos: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->SetYTitle("X Position (meters)");
  htemp->SetXTitle("Y Position (meters)");
  htemp->SetZTitle("RMS of Depth (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/RMS_plot_3d.jpg", imagesdirname);
  c2->SaveAs(filename);
  sprintf(filename, "%s/RMS_plot_3d.root", imagesdirname);
  c2->SaveAs(filename);

  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;
  
  return 0;
}