{

#include <iostream>
#include <fstream>
#include <string>
 
  gROOT->Reset();
  using namespace std;

  TCanvas *c1 = new TCanvas;
  gStyle->SetFrameFillStyle(0);
  gStyle->SetPadColor(0);
  gStyle->SetFrameFillColor(0);
  gStyle->SetCanvasColor(0);
  gStyle->SetOptStat(0);
  c1->SetTitle("");
  c1->SetLogx();
  c1->SetLogy();

  enum {NCOLORS = 240};
  int mycolors[NCOLORS];
  Float_t r, g, b;
  int i;
  for (i = 0; i<NCOLORS; i++){
    r = 1.0 - float(i)/float(NCOLORS);
    g = 1.0 - float(i)/float(NCOLORS);
    b = 1.0 - float(i)/float(NCOLORS);
    mycolors[i] = TColor::GetColor(r,g,b);
  }
  gStyle->SetPalette(NCOLORS, mycolors);
  
  double sfr_axis[500];
  
  double line;
  int loop = 0;
  
  ifstream sfr_axis_file("sfr_axis.txt");
  while (! sfr_axis_file.eof() ) {
    sfr_axis_file>>line;
    sfr_axis[loop] = line;
    loop++;
  }
  
  loop = 0;
  double age_axis[34];
  ifstream age_axis_file("model_age.txt");
  while (! age_axis_file.eof() ) {
    age_axis_file>>line;
    age_axis[loop] = line;
    loop++;
  }
  
  int loopx = 0;
  int loopy = 0;
  double sfrcomp[34][500];
  
  ifstream sfrcom_file("fpscomprobsfrhist.txt");
  while (! sfrcom_file.eof() ) {
    sfrcom_file>>line;
    sfrcomp[loopx][loopy] = line;
    //loop_x++;
    loopy++;
    if (loopy == 500){
      loopy = 0;
      loopx++;
    }
    if (loopx == 34)
      loopx = 0;
  }
  

  Double_t xmin = 1.e6;
  Double_t xmax = 2.8e9;
  Int_t xnbins = 20;
  Double_t logxmin = TMath::Log10(xmin);
  Double_t logxmax = TMath::Log10(xmax);
  Double_t xbinwidth = (logxmax - logxmin)/xnbins;
  Double_t xbins[xnbins+1];
  xbins[0] = xmin;

  for (i = 1; i<=xnbins; i++){
    xbins[i] =  TMath::Power(10,logxmin+i*xbinwidth);
  }


  Double_t ymin = 1.0;
  Double_t ymax = 6.e5;
  Int_t ynbins = 20;
  Double_t logymin = TMath::Log10(ymin);
  Double_t logymax = TMath::Log10(ymax);
  Double_t ybinwidth = (logymax - logymin)/ynbins;
  Double_t ybins[ynbins+1];
  ybins[0] = ymin;

  for (i = 1; i<=ynbins; i++){
    ybins[i] =  TMath::Power(10,logymin+i*ybinwidth);
  }




  TH2F *sfr_prob_hist = new TH2F("sfr_prob_hist","sfr_prob_hist",xnbins,xbins,ynbins,ybins);
  
  int j;
  float summerA =0.;
  float summer = 0.;
  double agehere = 0.;
  double sfrhere = 0.;
  double probhere = 0.;
  for (i = 0; i<34; i++) {
    for (j = 0; j<500; j++) {
      agehere = (age_axis[i]);
      sfrhere = (sfr_axis[j]);
      probhere = sfrcomp[i][j];
      summerA = summerA + probhere;
      sfr_prob_hist->Fill(agehere,sfrhere,probhere);
    }
  }
  //  sfr_prob_hist->GetXaxis()->SetBinLabel(1,"10^{6}");
  sfr_prob_hist->Draw("colz");
  
  // Int_t npeaks = 20;
  //  TSpectrum *s = new TSpectrum(2*npeaks);
  //Int_t nfound = s->Search(sfr_prob_hist,1,"");

  //make a 1D arrays of the sfr histogram
  int valnum = xnbins*ynbins;
  float age_holder[valnum];
  float sfr_holder[valnum];
  float prob_holder[valnum];
  int position_holder[valnum];

  int poscounter =0;
  float sumhere = 0;
  for (i=0; i<xnbins; i++){
    for (j = 0; j<ynbins; j++){
      age_holder[poscounter] = sfr_prob_hist->GetXaxis()->GetBinCenter(i);
      sfr_holder[poscounter] = sfr_prob_hist->GetYaxis()->GetBinCenter(j);
      prob_holder[poscounter] = sfr_prob_hist->GetBinContent(i,j);
      position_holder[poscounter] = poscounter;
      sumhere = sumhere + prob_holder[poscounter];
      // cout << sumhere << " " << poscounter << endl;
      poscounter++;
    }
  }
  
  //implement 1D bubble sort of prob_holder and keep track of bin position with position_holder

  Int_t flag = 1;
  float temp_hold;
  int temp_pos;
  int numLength = valnum;
  for (i=1;(i<valnum && flag); i++){
    flag = 0;
    for (j=0; j<(valnum-i); j++) {
      if (prob_holder[j+1]>prob_holder[j]){
	temp_hold = prob_holder[j];
	prob_holder[j] = prob_holder[j+1];
	prob_holder[j+1] = temp_hold;
	flag = 1;   // indicate a swap occurred

	//keep track of the position
	temp_pos = position_holder[j];
	position_holder[j] = position_holder[j+1];
	position_holder[j+1] = temp_pos;	
      }

    }
  }

  float summed_prob = 0.;
  int point_68 = 0;
  int flag_68 = 0;
  int point_95 = 0;
  int flag_95 = 0;
  // take sorted array and find the 68% and 95% positions
  for (i=0; i<valnum;i++) {
    summed_prob = summed_prob + prob_holder[i];
    //cout << summed_prob << " " << i << endl;
    if (summed_prob > 0.683 && flag_68==0) {
      point_68 = i;
      flag_68 = 1;
    }
    if (summed_prob > 0.954 && flag_95==0) {
      point_95 = i;
      flag_95 = 1;
    }
  }

  TH2F *sfr_prob_68 = new TH2F("sfr_prob_68","sfr_prob_68",xnbins,xbins,ynbins,ybins);

  agehere = 0.;
  sfrhere = 0.;
  probhere = 0.;
  int real_pos;
  poscounter = 0;
  for (i = 0; i<point_68; i++) {
    real_pos = position_holder[i];
    agehere = (age_holder[real_pos]);
    sfrhere = (sfr_holder[real_pos]);
    probhere = prob_holder[i];
    summer = summer + probhere;
    //cout << " real_pos is: " << real_pos << " probhere is: " << probhere << endl;
    sfr_prob_68->Fill(agehere,sfrhere,probhere);
    poscounter++;
  }

  //  sfr_prob_68->SetContour(2);
  //sfr_prob_68->SetContourLevel(0,0.01);
  //sfr_prob_68->Smooth(1,"k3a");
  c1->cd();
  sfr_prob_68->Draw("contz,list");
  c1->Update();

  TCanvas *c2 = new TCanvas("c2","First contour",100,100,800,600);
  TObjArray *contours = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
  if (!contours){
    cout << "entered here" << endl;
    return;
  }
  TList *lcontour1 = (TList*)contours->At(0);
  if (!lcontour1) return;
  TGraph *gcl1 = (TGraph*)lcontour1->First();
  if (!gcl1) return;
  if (gcl1->GetN() < 10) return;

  TList *lcontour2 = (TList*)contours->At(1);
  if (!lcontour2) return;
  TGraph *gcl2 = (TGraph*)lcontour2->First();
  if (!gcl2) return;
  if (gcl2->GetN() < 10) return;


  //gcl->SetMarkerStyle(21);
  //gcl->Draw("ac");
  int num_points = gcl1->GetN();
  Double_t tg_68_x[num_points];
  Double_t tg_68_y[num_points];
  float graph_68_x[num_points];
  float graph_68_y[num_points];
  for (i=0;i<num_points;i++){
    gcl1->GetPoint(i,tg_68_x[i],tg_68_y[i]);
    graph_68_x[i] = TMath::Power(10,tg_68_x[i]);
    graph_68_y[i] = TMath::Power(10,tg_68_y[i]);
  }

  TGraph *graph_68 = new TGraph(num_points,graph_68_x,graph_68_y);

}