void drawSlices(TH2F* h) {
  int nBinsX = h->GetNbinsX();
  int nBinsY = h->GetNbinsY();
  double rangeX[2] = {h->GetXaxis()->GetXmin(), h->GetXaxis()->GetXmax()};
  double rangeY[2] = {h->GetYaxis()->GetXmin(), h->GetYaxis()->GetXmax()};
  int binRangeX[2] = {h->GetXaxis()->GetFirst(), h->GetXaxis()->GetLast()};
  int binRangeY[2] = {h->GetYaxis()->GetFirst(), h->GetYaxis()->GetLast()};

//   std::vector<TH1F*> slicesX;
//   std::vector<TH1F*> slicesY;
  TCanvas* cY = new TCanvas("cY", "cYslices", 600, 550, 700, 400);
  cY->Divide(3,4, 0.0001, 0.0001);
  for (int iX = 1; iX <= nBinsX; ++iX) {
    TH1F* hslice = (TH1F*) h->ProjectionY(Form("%d_py", iX), iX, iX);
    //slicesX.push_back(hslice);
    cY->cd(iX);
    hslice->Draw();
    if (hslice->GetEntries()>0)
      printf("      slice #%d  N %5d, %10.4f +/- %10.4f\n", iX, hslice->GetEntries(), hslice->GetMean(), hslice->GetRMS());
  }
  for (int iY = 1; iY <= nBinsY; ++iY) {
  }
}

void printHistVals(TH1F* h) {
  if (h == NULL) return;
  int bins = h->GetNbinsX();
  for (int iBin = 1; iBin <= bins; ++iBin) {
    double content = h->GetBinContent(iBin);
    double error   = h->GetBinError(iBin);
    if (content != 0.0) {
      printf("      #%d  %10.4f +/- %10.4f\n", iBin, content, error);
    }
  }
}

void doProfiles(TH2F* h) {
  // Profile:
  // Calculates the Mean and RMS for each slice
  // Creates histograms where each bin is the mean value and its error is RMS/sqrt(slice entries)
  TCanvas* c = new TCanvas("cProfiles", "cProfiles", 0,0,1280, 1024);
  c->Divide(2,2,0.001, 0.001);
  c->cd(1);
  h->Draw("COLZ");
  c->cd(2);
  TH1F* hprofX = (TH1F*) h->ProfileY();
  hprofX->Draw();
  c->cd(3);
  TH1F* hprofY = (TH1F*) h->ProfileX();
  hprofY->Draw();

  printf("hprofX\n");
  printHistVals(hprofX);
  printf("hprofY\n");
  printHistVals(hprofY);
}

void doFitSlices(TH2F* h) {
  TCanvas* c = new TCanvas("cFitSlices", "cFitSlices", 0,0,1280, 1024);

  TF1* sliceFit[2];
  sliceFit[0] = new TF1("gausX", "gaus");
  sliceFit[1] = new TF1("gausY", "gaus");
  sliceFit[0]->SetParameters(100, 0.1, 0.5);
  sliceFit[1]->SetParameters(100, -0.5, 2.0);

  TObjArray hSlicesX;
  TObjArray hSlicesY;
  c->Divide(2,2,0.001, 0.001);
  c->cd(1);
  h->Draw("COLZ");
  c->cd(2);
  h->FitSlicesY(sliceFit[0], 0, -1, 0, "QNR", &hSlicesX);
  TH1F* fitSlicesX[4];
  int nFitPars = hSlicesX.GetEntries();
  for (int iPar = 0; iPar < nFitPars; ++iPar) {
    fitSlicesX[iPar] = (TH1F*) ((TH1F*) hSlicesX.At(iPar))->Clone();
    printf("Entries 1 %x  %d %s\n", hSlicesX.At(iPar), ((TH1F*) hSlicesX.At(iPar))->GetEntries(), hSlicesX.At(iPar)->GetName() );
    printf("Entries 2 %x  %d %s\n", fitSlicesX[iPar], fitSlicesX[iPar]->GetEntries(), fitSlicesX[iPar]->GetName() );
  }
  fitSlicesX[1]->Draw();

  c->cd(3);
//   TH1F* hsliceFitY = (TH1F*) h->FitSlicesY(sliceFit[1], 0, -1, 0, "QNR", NULL);//&aSlicesY);
//   TH1F* fitSlicesX[4];
//   fitSlicesY[0] = (TH1F*) gDirectory->Get("h_0")->Clone();
//   fitSlicesY[1] = (TH1F*) gDirectory->Get("h_1")->Clone();
//   fitSlicesY[2] = (TH1F*) gDirectory->Get("h_2")->Clone();
//   fitSlicesY[3] = (TH1F*) gDirectory->Get("h_chi2")->Clone();
//   fitSlicesY[1]->Draw();

  printf("hsliceFitX\n");
  printHistVals(fitSlicesX[1]);
  printf("hsliceFitY\n");
//  printHistVals(fitSlicesY[1]);
}

void profileHist(){
  TH1::SetDefaultSumw2();
  gStyle->SetPalette(1);
  TH2F* h = new TH2F("h", "ht", 10, -5, 5, 10, -10, 10);
  TF2 g2("gaus2d", "xgaus(0)*ygaus(3)");
  g2.SetParameters(100, 0.1, 0.5, 100, -0.5, 2);
  h->FillRandom("gaus2d", 10000);

  //doProfiles(h);
  doFitSlices(h);
  
  drawSlices(h);
}