#include "TH2F.h"
#include "TH1F.h"
#include "TF1.h"
#include "TROOT.h"
#include "TAxis.h"
#include "TStyle.h"
#include "TSystem.h"
#include "TInterpreter.h"
#include "TPad.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TMath.h"

TCanvas *c2;

Double_t myfunction1(Double_t *x, Double_t *par)   {
  Float_t xx =x[0];
  // Double_t f = par[3]+par[4]*xx+par[5]*xx*xx;
  Double_t f = par[4]*0.5*TMath::Exp(-xx*xx/(par[3]*par[3]));
  if (f<0) f=0;
  Double_t sigma=par[2];
  Double_t arg=(xx-par[1]);
  if (sigma==0) sigma=1e-19;
  arg/=sigma;
  f+=par[0]*TMath::Exp(-0.5*arg*arg);
  return f;
}

Double_t myfunction(Double_t *x, Double_t *par)   {
  Float_t xx =x[0];
  // Double_t f = par[3]+par[4]*xx+par[5]*xx*xx;
  Double_t f = par[0]*0.7*(1-(xx-90)*(xx-90)/(par[1]-90)/(par[1]-90));
  if (f<0) f=0;
  Double_t sigma=par[2];
  Double_t arg=(xx-par[1]);
  if (sigma==0) sigma=1e-19;
  arg/=sigma;
  f+=par[0]*TMath::Exp(-0.5*arg*arg);
  return f;
}
void myfunc()
{
  TF1 *f1 = new TF1("myfunc",myfunction,0,1000,3);
  f1->SetParameters(100,500,50,300,60);
  f1->SetParNames("constant","Mean","Sigma","p0","p1","p2");
  f1->Draw();
}
void myfit()
{
  TH1F *h1=new TH1F("h1","test",100,0,10);
  h1->FillRandom("myfunc",20000);
  TF1 *f1=(TF1 *)gROOT->GetFunction("myfunc");
  f1->SetParameters(800,1);
  h1->Fit("myfunc");
}


void MFitSlices(TH2 *h,int onX,TF1 *f1=0, Int_t firstbin=0, Int_t lastbin=0, Int_t cut=0, Option_t *option="", TObjArray* arr=NULL)
{
   TAxis& outerAxis = (onX ? *(h->GetYaxis()) : *(h->GetXaxis()));
   TAxis& innerAxis = (onX ? *(h->GetXaxis()) : *(h->GetYaxis()));

   Int_t nbins  = outerAxis.GetNbins();
   if (firstbin < 0) firstbin = 0;
   if (lastbin < 0 || lastbin > nbins + 1) lastbin = nbins + 1;
   if (lastbin < firstbin) {firstbin = 0; lastbin = nbins + 1;}
   TString opt = option;
   opt.ToLower();
   Int_t ngroup = 1;
   if (opt.Contains("g2")) {ngroup = 2; opt.ReplaceAll("g2","");}
   if (opt.Contains("g3")) {ngroup = 3; opt.ReplaceAll("g3","");}
   if (opt.Contains("g4")) {ngroup = 4; opt.ReplaceAll("g4","");}
   if (opt.Contains("g5")) {ngroup = 5; opt.ReplaceAll("g5","");}

   //default is to fit with a gaussian
   if (f1 == 0) {
      f1 = (TF1*)gROOT->GetFunction("gaus");
      if (f1 == 0) f1 = new TF1("gaus","gaus",innerAxis.GetXmin(),innerAxis.GetXmax());
      else         f1->SetRange(innerAxis.GetXmin(),innerAxis.GetXmax());
   }
   Int_t npar = f1->GetNpar();
   if (npar <= 0) return;
   Double_t *parsave = new Double_t[npar];
   f1->GetParameters(parsave);

   if (arr) {
      arr->SetOwner();
      arr->Expand(npar + 1);
   }

   //Create one histogram for each function parameter
   Int_t ipar;
   TH1D **hlist = new TH1D*[npar];
   char *name   = new char[2000];
   char *title  = new char[2000];
   const TArrayD *bins = outerAxis.GetXbins();
   for (ipar=0;ipar<npar;ipar++) {
      sprintf(name,"%s_%d",h->GetName(),ipar);
      sprintf(title,"Fitted value of par[%d]=%s",ipar,f1->GetParName(ipar));
      delete gDirectory->FindObject(name);
      if (bins->fN == 0) {
         hlist[ipar] = new TH1D(name,title, nbins, outerAxis.GetXmin(), outerAxis.GetXmax());
      } else {
         hlist[ipar] = new TH1D(name,title, nbins,bins->fArray);
      }
      hlist[ipar]->GetXaxis()->SetTitle(outerAxis.GetTitle());
      if (arr)
         (*arr)[ipar] = hlist[ipar];
   }
   sprintf(name,"%s_chi2",h->GetName());
   delete gDirectory->FindObject(name);
   TH1D *hchi2 = 0;
   if (bins->fN == 0) {
      hchi2 = new TH1D(name,"chisquare", nbins, outerAxis.GetXmin(), outerAxis.GetXmax());
   } else {
      hchi2 = new TH1D(name,"chisquare", nbins, bins->fArray);
   }
   hchi2->GetXaxis()->SetTitle(outerAxis.GetTitle());
   if (arr)
      (*arr)[npar] = hchi2;

   //Loop on all bins in Y, generate a projection along X
   Int_t bin;
   Long64_t nentries;
   TF1 *fgauss = new TF1("fgauss","gaus",0,100);

   for (bin=lastbin;bin>=firstbin;bin -= ngroup) {
      TH1D *hp;
      if (onX)
         hp= h->ProjectionX("_temp",bin,bin+ngroup-1,"");
      else
         hp= h->ProjectionY("_temp",bin,bin+ngroup-1,"");
      if (hp == 0) continue;
      nentries = Long64_t(hp->GetEntries());
      if (nentries == 0 || nentries < cut) {delete hp; continue;}
      TF1 *f0 = f1;
      //if (bin<5) f0=fg;
      //f0->SetParameters(parsave);
      

      hp->Fit(f0,opt.Data());
      // if (bin>2){
	fgauss->SetParameter(0,f0->GetParameter(0));
	fgauss->SetParameter(1,f0->GetParameter(1));
	fgauss->SetParameter(2,f0->GetParameter(2));
	
	fgauss->SetRange(f0->GetParameter(1)-5*f0->GetParameter(2)
	              ,f0->GetParameter(1)+3*f0->GetParameter(2));
	// }
      f0->SetRange(80,120+(outerAxis.GetBinCenter(bin)-3.5)*800/5);

      
      hp->GetYaxis()->SetRangeUser(0,f0->GetParameter(0)*2);
      //  hp->GetXaxis()->SetRangeUser(65,f0->GetParameter(1)*1.5);
      hp->GetXaxis()->SetRangeUser(60,120+(outerAxis.GetBinCenter(bin)-3.5)*800/5);
      if (bin<5){
	hp->GetYaxis()->SetRangeUser(0,2500/bin);
	//hp->GetXaxis()->SetRangeUser(80,200);
      }
      sprintf(name,"STM Bias %2.1f V;ADC;N",outerAxis.GetBinCenter(bin)+25);
      hp->SetTitle(name);
      sprintf(name,"%s_%2.1f.eps",h->GetName(),outerAxis.GetBinCenter(bin)+25);

      hp->Fit(fgauss,"RB+");
      hp->DrawCopy();

      c2->Update();
      //c2->SaveAs(name);
      getchar();
      Int_t npfits = f0->GetNumberFitPoints();
      //      if (npfits > npar && npfits >= cut) {
         Int_t binOn = bin + ngroup/2;
         for (ipar=0;ipar<npar;ipar++) {
            hlist[ipar]->Fill(outerAxis.GetBinCenter(binOn),fgauss->GetParameter(ipar));
            hlist[ipar]->SetBinError(binOn,fgauss->GetParError(ipar));
         }
         hchi2->Fill(outerAxis.GetBinCenter(binOn),f0->GetChisquare()/(npfits-npar));
	 printf("%2.1f %2.1f  %2.1f %2.1f\n",outerAxis.GetBinCenter(binOn),fgauss->GetParameter(0),fgauss->GetParameter(1),fgauss->GetParameter(2));

	 //}
	 //      delete hp;
   }
   delete [] parsave;
   delete [] name;
   delete [] title;
   delete [] hlist;
}

void fitslicesy() {
//
// To see the output of this macro, click begin_html <a href="gif/fitslicesy.gif" >here</a> end_html
//    This macro illustrates how to use the TH1::FitSlicesY function
//    It uses the TH2F histogram generated in macro hsimple.C
//    It invokes FitSlicesY and draw the fitted "mean" and "sigma"
//    in 2 sepate pads.
//    This macro shows also how to annotate a picture, change
//    some pad parameters.
//Author: Rene Brun
   
// Change some default parameters in the current style
   gStyle->SetLabelSize(0.06,"x");
   gStyle->SetLabelSize(0.06,"y");
   gStyle->SetFrameFillColor(kWhite);
   gStyle->SetTitleW(0.6);
   gStyle->SetTitleH(0.05);
   gStyle->SetOptFit(0);
   gStyle->SetOptStat(0);

// Connect the input file and get the 2-d histogram in memory
   TString dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName());
   dir.ReplaceAll("fitslicesy.C","hsimple.C");
   dir.ReplaceAll("/./","/");
   if (!gInterpreter->IsLoaded(dir.Data())) gInterpreter->LoadMacro(dir.Data());
   TFile *hsimple =  new TFile("stm_bias_scan0.root");
     //TFile *hsimple = (TFile*)gROOT->ProcessLineFast("hsimple(1)");
   if (!hsimple) return;
   hsimple->ls();
   //TH2F *hpxpy = (TH2F*)hsimple->Get("hpxpy");
   TH2F *hpxpy = (TH2F*)hsimple->Get("h2dadcdac001");
   if (!hpxpy) return;
   c2 = new TCanvas("c2","c2",700,500);
   c2->cd();
// Fit slices projected along Y fron bins in X [7,32]
   //TF1 *fb = new TF1("fb","gaus",140,1000);
   TF1 *fb = new TF1("myfunc",myfunction,90,1000,3);
   fb->SetParameters(80,700,20);
   fb->SetParLimits(0,50,9000);
   fb->SetParLimits(1,100,800);
   fb->SetParLimits(2,10,80);
   fb->SetParNames("constant","Mean","Sigma","p0","p1","p2");

   TF1 *fc = new TF1("fc","gaus",0,110);
   // fitsliy(hpxpy,fb);
 
   MFitSlices(hpxpy,1,fc,1,1000,-1120,"RBQ");
   c2->Update();
   

// Create a canvas and divide it
   TCanvas *c1 = new TCanvas("c1","c1",700,500);
   c1->SetFillColor(42);
   c1->Divide(2,1);
   c1->cd(1);
   TPad *left = (TPad*)gPad;
   left->Divide(1,2);

// Draw 2-d original histogram
   left->cd(1);
   gPad->SetTopMargin(0.12);
   gPad->SetFillColor(33);
   hpxpy->Draw("col");
   hpxpy->GetXaxis()->SetLabelSize(0.06);
   hpxpy->GetYaxis()->SetLabelSize(0.06);
   hpxpy->SetMarkerColor(kYellow);


   char name[128];
  
 sprintf(name,"%s_chi2",hpxpy->GetName());  TH1D *hpxpy_0=gDirectory->FindObject(name);
 sprintf(name,"%s_1",hpxpy->GetName());  TH1D *hpxpy_1=gDirectory->FindObject(name);
 sprintf(name,"%s_2",hpxpy->GetName());  TH1D *hpxpy_2=gDirectory->FindObject(name);
// Show fitted "mean" for each slice



   left->cd(2);
   gPad->SetFillColor(33);
   hpxpy_0->Draw();
   c1->cd(2);
   TPad *right = (TPad*)gPad;
   right->Divide(1,2);
   right->cd(1);
   gPad->SetTopMargin(0.12);
   gPad->SetLeftMargin(0.15);
   gPad->SetFillColor(33);
   hpxpy_1->Draw();

// Show fitted "sigma" for each slice
   right->cd(2);
   gPad->SetTopMargin(0.12);
   gPad->SetLeftMargin(0.15);
   gPad->SetFillColor(33);
   hpxpy_2->SetMinimum(0.8);
   hpxpy_2->Draw();

//attributes
   hpxpy_0->SetLineColor(kYellow);
   hpxpy_1->SetLineColor(kYellow);
   hpxpy_2->SetLineColor(kYellow);
   hpxpy_0->SetMarkerColor(kRed);
   hpxpy_1->SetMarkerColor(kRed);
   hpxpy_2->SetMarkerColor(kRed);
   hpxpy_0->SetMarkerStyle(21);
   hpxpy_1->SetMarkerStyle(21);
   hpxpy_2->SetMarkerStyle(21);
   hpxpy_0->SetMarkerSize(0.6);
   hpxpy_1->SetMarkerSize(0.6);
   hpxpy_2->SetMarkerSize(0.6);
}