#include <iostream>

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

char rootfilename[] = "zygo.root";

void t(char *upper_path="E6-10_2take2/stch_upper")
{
   TFile *rfile = new TFile(rootfilename);
   if (! rfile->IsOpen()) {
      std::cerr << "Error - cannot open input root file, giving up"
                << std::endl;
      return;
   }
   TH2D *hsu;
   Map2D *su=0;
   TDirectory *du = rfile->GetDirectory(upper_path,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");
   }

   su->SetContour(80);
   TCanvas *c1 = (TCanvas*)gROOT->FindObject("c1");
   if (c1 == 0) {
      c1 = new TCanvas("c1","c1",500,500);
   }
   else {
      c1->cd();
   }
   su->Draw("surf4");
   TCanvas *c2 = (TCanvas*)gROOT->FindObject("c2");
   if (c2 == 0) {
      c2 = new TCanvas("c2","c2",500,500);
   }
   else {
      c2->cd();
   }
   gStyle->SetPalette(1);
   su->Draw("colz");
}

void t2(char *upper_path="E6-10_2take2/stch_upper",
        char *lower_path="E6-10_2take2/stch_lower")
{
   TFile *rfile = new TFile(rootfilename);
   if (! rfile->IsOpen()) {
      std::cerr << "Error - cannot open input root file, giving up"
                << std::endl;
      return;
   }
   TH2D *hsu;
   Map2D *su=0;
   TDirectory *du = rfile->GetDirectory(upper_path,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;
   TDirectory *dl = rfile->GetDirectory(lower_path,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");
   }

   sl->Shift(0,0,+40);

   double refindex=2.42;
   Map2D *rdiff = new Map2D(*oru);
   rdiff->Add(orl,-1);
   rdiff->Shift(0,0,-40);
   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 *newlow = new Map2D(*su);
   newlow->Add(thick,-1);
   newlow->SetContour(80);
   newlow->SetName("newlow");
   newlow->SetTitle("lower surface, corrected");  
   TCanvas *c1 = (TCanvas*)gROOT->FindObject("c1");
   if (c1 == 0) {
      c1 = new TCanvas("c1","c1",500,500);
   }
   else {
      c1->cd();
   }
   newlow->Draw("surf4");
   TCanvas *c2 = (TCanvas*)gROOT->FindObject("c2");
   if (c2 == 0) {
      c2 = new TCanvas("c2","c2",500,500);
   }
   else {
      c2->cd();
   }
   gStyle->SetPalette(1);
   thick->Draw("colz");
}