//Example for fitting signal/background. 
// This example can be executed with:
// root > .x FittingDemo.C  (using the CINT interpreter)
// root > .x FittingDemo.C+ (using the native complier via ACLIC)
//Author: Rene Brun
   
#include "TH1.h"
#include "TMath.h"
#include "TF1.h"
#include "TLegend.h"
#include "TCanvas.h"

// Quadratic background function
Double_t background(Double_t *x, Double_t *par) {
   return par[0] + par[1]*x[0] + par[2]*x[0]*x[0];
}


// Lorenzian Peak function
Double_t lorentzianPeak(Double_t *x, Double_t *par) {
  return (0.5*par[0]*par[1]/TMath::Pi()) / 
    TMath::Max( 1.e-10,(x[0]-par[2])*(x[0]-par[2]) 
   + .25*par[1]*par[1]);
}

// Sum of background and peak function
Double_t fitFunction(Double_t *x, Double_t *par) {
  return background(x,par) + lorentzianPeak(x,&par[3]);
}

void NewFittingDemo() {
 //Bevington Exercise by Peter Malzacher, modified by Rene Brun
 
   ifstream in;
   in.open("./output.txt");
   Int_t x,y, data[1000000][2], Hash[1000000];
   Double_t xAry[1000000], yAry[1000000];
   Long64_t index[1000000], n=0;
   Double_t yMean=0,xMean=0, XShift;
   //Choose the mean 0:the first mean, 1: the second and so on
   int Mean=1; 
   for(Int_t i=0; i<1000000; i++)
      Hash[i]=0;
   while (1) {
     in >> x;
     if (!in.good()) break;
     Hash[x]++;
    }
    
   for(i=0; i<1000000; i++)
    {
      if (Hash[i]!=0)
       {
        data[n][0]=i-150;
        data[n][1]=Hash[i]; 
        n++;
       }
    }
  
   in.close();
   Int_t p1=0;
   //the peak on intervals- set the size of interval on x axis
   XShift=1000;
   for(Int_t i=0; (i+XShift)<n; i+=XShift)
   {
   	   yMean=0;
   	   for (r=0; r<n; r++)
   	   {
   	      if((data[r][0]>=i)&&(data[r][0]<=i+XShift))
   	       {
   	         if(data[r][1]>yMean))
   	          {
   	      	   yMean=data[r][1];
   	      	   xMean=data[r][0];
   	          }
   	          
   	       }
   	    }
   	   index[p1]=p1;
   	   yAry[p1]=yMean;
   	   xAry[p1]=xMean;
   	   p1++;
   }
   TMath::Sort(p1, yAry, index, 1);
   const int nBins = n;


   
  TCanvas *c1 = new TCanvas("c1","Fitting Demo",10,10,700,500);
   
   //c1->SetFillColor(30);
   //c1->SetFrameFillColor(41);
   c1->SetGrid();
   
   TH1F *histo = new TH1F("histo","
Peak Fitting",8192,0,8192);
   histo->SetMarkerStyle(20);
   histo->SetMarkerSize(0.8);
   histo->SetStats(0);
       
   for(int i=0; i < nBins;  i++) histo->SetBinContent(data[i][0],data[i][1]);
   
   // create a TF1 with the range from 0 to 3 and 6 parameters
   TF1 *fitFcn = new TF1("fitFcn",fitFunction,0,8192,6);
   fitFcn->SetNpx(500);
   fitFcn->SetLineWidth(4);
   fitFcn->SetLineColor(kMagenta);
    
   // first try without starting values for the parameters
   // This defaults to 1 for each param. 
   // this results in an ok fit for the polynomial function
   // however the non-linear part (lorenzian) does not 
   // respond well.
   //fitFcn->SetParameters(1,1,1,1,1,1);
   histo->Fit("fitFcn","0");
   
   // second try: set start values for some parameters
   fitFcn->SetParameter(4,80); // width
   fitFcn->SetParameter(5,xAry[index[Mean]]);   // peak
 
   histo->Fit("fitFcn","V+","ep");
   
   // improve the picture:
   TF1 *backFcn = new TF1("backFcn",background,0,8192,3);
   backFcn->SetLineColor(kRed);
   TF1 *signalFcn = new TF1("signalFcn",lorentzianPeak,0,8192,3);
   signalFcn->SetLineColor(kBlue);
   signalFcn->SetNpx(500);
   Double_t par[6];
  
   // writes the fit results into the par array
   fitFcn->GetParameters(par);
   backFcn->SetParameters(par);
   backFcn->Draw("same");
 
   signalFcn->SetParameters(&par[3]);
   signalFcn->Draw("same"); 
   
   // draw the legend
   TLegend *legend=new TLegend(0.6,0.65,0.88,0.85);
   legend->SetTextFont(72);
   legend->SetTextSize(0.04);
   legend->AddEntry(histo,"Data","lpe");
   legend->AddEntry(backFcn,"Background fit","l");
   legend->AddEntry(signalFcn,"Signal fit","l");
   legend->AddEntry(fitFcn,"Global Fit","l");
   legend->Draw();
   printf("=============display peaks==============\n");
   for(i=0; i<p1; i++)
     printf("rank=%d   xMean=%f    Peak=%f\n", i, xAry[index[i]], yAry[index[i]]);   
   printf("--------------------------------------------\n");
   printf("-----------------Mean=%f----------------\n",xAry[index[Mean]]);
   printf("-----------------Number of Peaks=%d---------------------------\n",p1);
   printf("-----------------Integral of Signal Function=%f-------------------\n",signalFcn->Integral(0,n));
 cout<<"\Integration="<<par[3]<<"  Mean="<<par[5]<<"  Sigma="<<par[4];  
}