// 
// plots.C - root macro to read Zygo image maps from a root file
// t
//           (previously created by zygo2root) and generate
//           pretty pictures of the upper and lower surfaces,
//           and the sample thickness profile, both with and
//           without stitching corrections.
//
// author: Richard Jones and Fridah Mokaya
// date: Sept. 26, 2012
//
// Usage example:
//	root [0] .L MetroProMap_cc.c+O
//	root [1] .L Map2D_cc.c+O
//      root [2] .x plots.C("myzygoplots.root","take1");
//      root [3] plotThickness();
//      root [5] plotUpperSurface();
//      root [7] plotLowerSurface();

double refindex=2.42;

#include <iostream>
#include <string>

#include <TROOT.h>
#include <TCanvas.h>
#include <TH2D.h>
#include "Map2D.h"

void plots(char *rootfile, char *rootpath)
{
   TFile *rfile = new TFile(rootfile);
   if (! rfile->IsOpen()) {
      std::cerr << "Error - cannot open input root file, giving up"
                << std::endl;
      return;
   }
   
   TH2D *hsu;
   Map2D *su=0;
   std::string upper_path(rootpath);
   upper_path += "/stch_upper";
   TDirectory *du = rfile->GetDirectory(upper_path.c_str(),kTRUE);
   du->GetObject("over2",hsu);
   if (hsu) {
      su = new Map2D(*hsu);
      su->SetName("upper");
      su->SetTitle("upper surface, stitched");
   }
   TH2D *hosu;
   Map2D *osu=0;
   du->GetObject("stitch_upper.dat",hosu);
   if (hosu) {
      osu = new Map2D(*hosu);
      osu->SetName("oupper");
      osu->SetTitle("upper surface, Zygo stitched");
   }
   TH2D *horu;
   Map2D *oru=0;
   du->GetObject("over1",horu);
   if (horu) {
      oru = new Map2D(*horu);
      oru->SetName("rupper");
      oru->SetTitle("upper surface, unstitched");
   }

   TH2D *hsl;
   Map2D *sl=0;
   std::string lower_path(rootpath);
   lower_path += "/stch_lower";
   TDirectory *dl = rfile->GetDirectory(lower_path.c_str(),kTRUE);
   dl->GetObject("over2",hsl);
   if (hsl) {
      sl = new Map2D(*hsl);
      sl->SetName("lower");
      sl->SetTitle("lower surface, stitched");
   }
   TH2D *hosl;
   Map2D *osl=0;
   dl->GetObject("stitch_lower.dat",hosl);
   if (hosl == 0) {
      dl->GetObject("stitch_lower2.dat",hosl);
   }
   if (hosl == 0) {
      dl->GetObject("stitch_lower3.dat",hosl);
   }
   if (hosl) {
      osl = new Map2D(*hosl);
      osl->SetName("olower");
      osl->SetTitle("lower surface, Zygo stitched");
   }
   TH2D *horl;
   Map2D *orl=0;
   dl->GetObject("over1",horl);
   if (horl) {
      orl = new Map2D(*horl);
      orl->SetName("rlower");
      orl->SetTitle("lower surface, unstitched");
   }

   Map2D *rdiff = new Map2D(*oru);
   rdiff->Add(orl,-1);
   rdiff->Rescale(1,1,1/refindex);
   rdiff->SetName("rdiff");
   rdiff->SetTitle("thickness, unstitched");

   Map2D *thick = new Map2D(*su);
   thick->Add(sl,-1);
   thick->Rescale(1,1,1/refindex);
   thick->SetName("thick");
   thick->SetTitle("thickness, stitched");

   Map2D *lowcor = new Map2D(*su);
   lowcor->Add(thick,-1);
   lowcor->SetName("lowcor");
   lowcor->SetTitle("corrected lower surface, stitched");  

   Map2D *rlowcor = new Map2D(*oru);
   rlowcor->Add(rdiff,-1);
   rlowcor->SetName("rlowcor");
   rlowcor->SetTitle("corrected lower surface, unstitched");  
}

void plotThickness(char *option="colz")
{
   TCanvas *c1 = (TCanvas*)gROOT->FindObject("c1");
   if (c1 == 0) {
      c1 = new TCanvas("c1","c1",500,500);
   }
   else {
      c1->cd();
   }
   Map2D *rdiff = gROOT->FindObject("rdiff");
   std::string options(option);
   if (options == "colz" || options == "COLZ") {
      c1->SetRightMargin(0.2);
      rdiff->SetContour(150);
      gStyle->SetPalette(1);
      rdiff->Draw("colz");
      c1->Update();
      TPaletteAxis *palette = (TPaletteAxis*)rdiff->GetListOfFunctions()->FindObject("palette");
      rdiff->GetXaxis()->SetTitle("x (mm)");
      rdiff->GetYaxis()->SetTitle("y (mm)");
      rdiff->GetXaxis()->SetLimits(0,5);
      rdiff->GetYaxis()->SetLimits(0,5);
      rdiff->GetYaxis()->SetTitleOffset(1.2);
      palette->GetAxis()->SetTitle("microns");
      palette->SetTitleOffset(1.2);
      palette->SetX1NDC(0.82);
      palette->SetX2NDC(0.86);
   }
   rdiff->Draw(option);

   TCanvas *c2 = (TCanvas*)gROOT->FindObject("c2");
   if (c2 == 0) {
      c2 = new TCanvas("c2","c2",500,500);
   }
   else {
      c2->cd();
   }
   Map2D *thick = gROOT->FindObject("thick");
   std::string options(option);
   if (options == "colz" || options == "COLZ") {
      c2->SetRightMargin(0.2);
      thick->SetContour(150);
      gStyle->SetPalette(1);
      thick->Draw("colz");
      c2->Update();
      TPaletteAxis *palette = (TPaletteAxis*)thick->GetListOfFunctions()->FindObject("palette");
      thick->GetXaxis()->SetTitle("x (mm)");
      thick->GetYaxis()->SetTitle("y (mm)");
      thick->GetXaxis()->SetLimits(0,5);
      thick->GetYaxis()->SetLimits(0,5);
      thick->GetYaxis()->SetTitleOffset(1.2);
      palette->GetAxis()->SetTitle("microns");
      palette->SetTitleOffset(1.2);
      palette->SetX1NDC(0.82);
      palette->SetX2NDC(0.86);
   }
   thick->Draw(option);
}

void plotUpperSurface(char *option="colz")
{
   TCanvas *c1 = (TCanvas*)gROOT->FindObject("c1");
   if (c1 == 0) {
      c1 = new TCanvas("c1","c1",500,500);
   }
   else {
      c1->cd();
   }
   Map2D *rupper = gROOT->FindObject("rupper");
   double midheight = rupper->GetMaximum();
   rupper->Shift(0,0,-(midheight-50.));
   std::string options(option);
   if (options == "colz" || options == "COLZ") {
      c1->SetRightMargin(0.2);
      rupper->SetContour(150);
      gStyle->SetPalette(1);
      rupper->Draw("colz");
      c1->Update();
      TPaletteAxis *palette = (TPaletteAxis*)rupper->GetListOfFunctions()->FindObject("palette");
      rupper->GetXaxis()->SetTitle("x (mm)");
      rupper->GetYaxis()->SetTitle("y (mm)");
      rupper->GetXaxis()->SetLimits(0,5);
      rupper->GetYaxis()->SetLimits(0,5);
      rupper->GetYaxis()->SetTitleOffset(1.2);
      palette->GetAxis()->SetTitle("microns");
      palette->SetTitleOffset(1.2);
      palette->SetX1NDC(0.82);
      palette->SetX2NDC(0.86);
   }
   rupper->Draw(option);

   TCanvas *c2 = (TCanvas*)gROOT->FindObject("c2");
   if (c2 == 0) {
      c2 = new TCanvas("c2","c2",500,500);
   }
   else {
      c2->cd();
   }
   Map2D *upper = gROOT->FindObject("upper");
   midheight = upper->GetMaximum();
   upper->Shift(0,0,-(midheight-50.));
   std::string options(option);
   if (options == "colz" || options == "COLZ") {
      c2->SetRightMargin(0.2);
      upper->SetContour(150);
      gStyle->SetPalette(1);
      upper->Draw("colz");
      c2->Update();
      TPaletteAxis *palette = (TPaletteAxis*)upper->GetListOfFunctions()->FindObject("palette");
      upper->GetXaxis()->SetTitle("x (mm)");
      upper->GetYaxis()->SetTitle("y (mm)");
      upper->GetXaxis()->SetLimits(0,5);
      upper->GetYaxis()->SetLimits(0,5);
      upper->GetYaxis()->SetTitleOffset(1.2);
      palette->GetAxis()->SetTitle("microns");
      palette->SetTitleOffset(1.2);
      palette->SetX1NDC(0.82);
      palette->SetX2NDC(0.86);
   }
   upper->Draw(option);
}

void plotLowerSurface(char *option="colz")
{
   TCanvas *c1 = (TCanvas*)gROOT->FindObject("c1");
   if (c1 == 0) {
      c1 = new TCanvas("c1","c1",500,500);
   }
   else {
      c1->cd();
   }
   Map2D *rlowcor = gROOT->FindObject("rlowcor");
   double midheight = rlowcor->GetMaximum();
   rlowcor->Shift(0,0,-(midheight-50));
   std::string options(option);
   if (options == "colz" || options == "COLZ") {
      c1->SetRightMargin(0.2);
      rlowcor->SetContour(150);
      gStyle->SetPalette(1);
      rlowcor->Draw("colz");
      c1->Update();
      TPaletteAxis *palette = (TPaletteAxis*)rlowcor->GetListOfFunctions()->FindObject("palette");
      rlowcor->GetXaxis()->SetTitle("x (mm)");
      rlowcor->GetYaxis()->SetTitle("y (mm)");
      rlowcor->GetXaxis()->SetLimits(0,5);
      rlowcor->GetYaxis()->SetLimits(0,5);
      rlowcor->GetYaxis()->SetTitleOffset(1.2);
      palette->GetAxis()->SetTitle("microns");
      palette->SetTitleOffset(1.2);
      palette->SetX1NDC(0.82);
      palette->SetX2NDC(0.86);
   }
   rlowcor->Draw(option);

   TCanvas *c2 = (TCanvas*)gROOT->FindObject("c2");
   if (c2 == 0) {
      c2 = new TCanvas("c2","c2",500,500);
   }
   else {
      c2->cd();
   }
   Map2D *lowcor = gROOT->FindObject("lowcor");
   midheight = lowcor->GetMaximum();
   lowcor->Shift(0,0,-(midheight-50));
   std::string options(option);
   if (options == "colz" || options == "COLZ") {
      c2->SetRightMargin(0.2);
      lowcor->SetContour(150);
      gStyle->SetPalette(1);
      lowcor->Draw("colz");
      c2->Update();
      TPaletteAxis *palette = (TPaletteAxis*)lowcor->GetListOfFunctions()->FindObject("palette");
      lowcor->GetXaxis()->SetTitle("x (mm)");
      lowcor->GetYaxis()->SetTitle("y (mm)");
      lowcor->GetXaxis()->SetLimits(0,5);
      lowcor->GetYaxis()->SetLimits(0,5);
      lowcor->GetYaxis()->SetTitleOffset(1.2);
      palette->GetAxis()->SetTitle("microns");
      palette->SetTitleOffset(1.2);
      palette->SetX1NDC(0.82);
      palette->SetX2NDC(0.86);
   }
   lowcor->Draw(option);
}