#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]) + lorentzianPeak(x,&par[6]);
}

void FittingDemo() {
 //Bevington Exercise by Peter Malzacher, modified by Rene Brun

   const int nBins = 60;
Double_t data[nBins] = { 6, 1, 10, 12, 6, 13, 23, 22, 15, 21, 23, 26, 36, 25, 27, 35, 40, 44, 
66, 81, 75, 57, 48, 45, 46, 41, 35, 36, 35, 32, 40, 50, 55, 65, 70, 
65, 55, 37, 32, 32, 32, 33, 35, 33, 33, 39, 29, 31, 32, 34, 26, 39, 
29, 35, 32, 21, 21, 15, 25, 15};

   TCanvas *c1 = new TCanvas("c1","Fitting Demo",10,10,700,500);
   c1->SetFillColor(33);
   c1->SetFrameFillColor(41);
   c1->SetGrid();

   TH1F *histo = new TH1F("histo",
      "Lorentzian Peak on Quadratic Background",60,0,3);
   histo->SetMarkerStyle(21);
   histo->SetMarkerSize(0.8);
   histo->SetStats(0);

   for(int i=0; i < nBins;  i++) histo->SetBinContent(i+1,data[i]);

   // create a TF1 with the range from 0 to 3 and 6 parameters
   TF1 *fitFcn = new TF1("fitFcn",fitFunction,0,3,9);
   fitFcn->SetNpx(500);
   fitFcn->SetLineWidth(4);
   fitFcn->SetLineColor(kMagenta);
   
   // set reasonable initial values for all parameters
   fitFcn->SetParameters(0., 40., -10., // backFcn
                         15., 0.2, 1., // signalFcn1
                         10, 0.15, 1.7); // signalFcn2
   histo->Fit("fitFcn","0");
   
   histo->Fit("fitFcn","V+","ep");
   
   // improve the picture:
   TF1 *backFcn = new TF1("backFcn",background,0,3,3);
   backFcn->SetLineColor(kRed);
   TF1 *signalFcn1 = new TF1("signalFcn1",lorentzianPeak,0,3,3);
   signalFcn1->SetLineColor(kBlue);
   signalFcn1->SetNpx(500);
   TF1 *signalFcn2 = new TF1("signalFcn2",lorentzianPeak,1.2,3,3);
   signalFcn2->SetLineColor(kBlack);
   signalFcn2->SetNpx(500);
 
   Double_t par[9];

   // writes the fit results into the par array
   fitFcn->GetParameters(par);

   backFcn->SetParameters(par);
   backFcn->Draw("same");

   signalFcn1->SetParameters(&par[3]);
   signalFcn1->Draw("same");
   
   signalFcn2->SetParameters(&par[6]);
   signalFcn2->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(signalFcn1,"Signal fit1","l");
   legend->AddEntry(signalFcn2,"Signal fit2","l");
   legend->AddEntry(fitFcn,"Global Fit","l");
   legend->Draw();

}