#include "/usr/local/include/headfile_loading.h"

//int main() {
int plot_topo_contour() {

  gROOT->SetStyle("Plain");
  gStyle->TStyle::SetNumberContours(100);
  TFile * fin = new TFile("derived_length.root");
  
  TH2D * h2d = (TH2D*)fin->Get("h2d_length");

  //h2d->Smooth();
  //h2d->Smooth();
  //h2d->Smooth();
  int Ntheta = h2d->GetNbinsY();
  int Nphi   = h2d->GetNbinsX();

  FILE * file;
  file = fopen("derived_length_C1.dat", "w");
  
  double * x = new double[Ntheta*Nphi];
  double * y = new double[Ntheta*Nphi];
  double * z = new double[Ntheta*Nphi];

  //double x_det = 0;
  //double y_det = 0;
  //double z_det = 0;
  
  double x_det = 5288.1190;
  double y_det = 3368.8013;
  double z_det = 3629.8090;
  
  int index = 0;
  for (int itheta=1; itheta<=Ntheta; itheta++) {

    if ( itheta*TMath::Pi()/2./Ntheta > 50./180*TMath::Pi() ) continue;
    
    for (int iphi=1; iphi<=Nphi; iphi++) {
      
      double l = h2d->GetBinContent(iphi, itheta);
      
      double theta = h2d->GetYaxis()->GetBinCenter(itheta);
      double phi   = h2d->GetXaxis()->GetBinCenter(iphi);
      
      x[index] = l * TMath::Sin(theta) * TMath::Cos(phi) + x_det;
      y[index] = l * TMath::Sin(theta) * TMath::Sin(phi) + y_det;
      z[index] = l * cos(theta) + z_det;
      
      if (fabs(x[index])>1e-10) {

	cout << "l: " <<  l << " theta=" << theta << " phi=" << phi << endl;
      
	// 	printf("%15f %15f %15f\n",
	// 	       x[index] + x_det,
	// 	       y[index] + y_det,
	// 	       z[index] + z_det);
	
	fprintf(file, "%f %f %f 0.0\n",
		x[index],
		y[index],
		z[index]);
	index++;

      }
    }
  }

  TCanvas * can = new TCanvas("can", "", 800, 600);
  can->SetRightMargin(0.15);
  
  TGraph2D * g2d = new TGraph2D(index, x, y, z);

  g2d->Draw("contz");
  g2d->Draw("cont3same");
  
  //g2d->Draw("surf7z");
  //g2d->Draw("cont1z");
  //g2d->Draw("cont3same list");
  //g2d->Draw("CONT1 same");

  g2d->GetZaxis()->SetRangeUser(z_det, 500+z_det);

  //TObjArray *contours = gROOT->GetListOfSpecials()->FindObject("contours");
  //Int_t ncontours     = contours->GetSize();
  //TList *list         = (TList*)contours->At(i);
  
  double xmin = g2d->GetXaxis()->GetXmin();
  double xmax = g2d->GetXaxis()->GetXmax();

  double ymin = g2d->GetYaxis()->GetXmin();
  double ymax = g2d->GetYaxis()->GetXmax();

  double zmin = g2d->GetZaxis()->GetXmin();
  double zmax = g2d->GetZaxis()->GetXmax();
  
  cout << endl;
  cout << xmin << "  " << xmax << endl;
  cout << ymin << "  " << ymax << endl;


  double * xx_det = new double[1];
  double * yy_det = new double[1];
  double * zz_det = new double[1];

  xx_det[0] = x_det;
  yy_det[0] = y_det;
  zz_det[0] = zmax;
  
  TGraph2D * g_det = new TGraph2D(1, xx_det, yy_det, zz_det);
  g_det->SetMarkerStyle(20);
  g_det->SetMarkerColor(kRed);
  //g_det->Draw("psame");
  
  can->Update();
  can->Modified();
  can->Print("derived_length_contour.png");

  
  
}