plot_topo_contour() {

   TFile * fin = new TFile("derived_length.root");
   TH2D * h2d = (TH2D*)fin->Get("h2d_length");

   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 = 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) {
            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);
   gStyle->SetPalette(1);
  
   TGraph2D * g2d = new TGraph2D(index, x, y, z);
   TH2D *h2 = g2d->GetHistogram();
   h2->SetMinimum(z_det);
   h2->SetMaximum(500+z_det);
  
   h2->Draw("cont4z");
   gPad->Update();

   Double_t bm = gPad->GetBottomMargin();
   Double_t lm = gPad->GetLeftMargin();
   Double_t rm = gPad->GetRightMargin();
   Double_t tm = gPad->GetTopMargin();
   Double_t x1 = h2->GetXaxis()->GetXmin();
   Double_t y1 = h2->GetYaxis()->GetXmin();
   Double_t x2 = h2->GetXaxis()->GetXmax();
   Double_t y2 = h2->GetYaxis()->GetXmax();

   TPad *padtrans = new TPad("padtrans","padtrans",0,0,1,1);
   padtrans->SetFillStyle(0);
   padtrans->SetFrameFillStyle(0);
   padtrans->Draw();
   padtrans->cd();
   padtrans->Range(x1-(x2-x1)*(lm/(1-rm-lm)),
                   y1-(y2-y1)*(bm/(1-tm-lm)),
		   x2+(x2-x1)*(rm/(1-rm-lm)), 
		   y2+(y2-y1)*(tm/(1-tm-lm)));
   gPad->Update();
   h2->DrawClone("cont3 same");
}