const unsigned int n_bins(100);
//const double m_min(0.1399), m_max(0.16);
const double m_min(0.1399), m_max(0.159);
const double bin_width = (m_max-m_min)/((double)n_bins);

double gaus_func(double m, double mu, double sigma)
{
	
  if (m<m_min || m>m_max) return 0.;
  double norm = TMath::Sqrt(TMath::Pi()/2.)*sigma*(TMath::Erf((mu-m_min)/(TMath::Sqrt2()*sigma))-TMath::Erf((mu-m_max)/(TMath::Sqrt2()*sigma)));
  if (norm==0.) return 0.;
  return TMath::Gaus(m,mu,sigma,0)/norm;
}

double sgn_func(double *x, double *p)
{
  double m = x[0];
  if (m<m_min || m>m_max) return 0.;
  
  double Nsig = p[0];
  double mu = p[1];
  double sigma= p[2];
    
  return Nsig*(gaus_func(m,mu,sigma))*bin_width;
}


double bkg_func(double *x, double *p)
{
  double m = x[0];
  if (m<m_min || m>m_max) return 0.;
  
  double Nbkg = p[0];
  double alpha = p[1];
  double beta = p[2];
  double m0 = p[3];
	
  
  double norm = ((2./3.) * (pow((m_max-m0),3./2.)) + (2./5.) * alpha* (pow((m_max-m0),5./2.)) + (2./7.) * beta *(pow((m_max-m0),7./2.))) - ((2./3.) * (pow((m_min-m0),3./2.)) +  (2./5.) * alpha* (pow((m_min-m0),5./2.)) + (2./7.) * beta *(pow((m_min-m0),7./2.)));
  
  
  if (norm==0.) return 0.;
  return Nbkg*((pow((m-m0),1./2.) + alpha*pow((m-m0),3./2.) + beta*pow((m-m0),5./2.))/norm)*bin_width;
  
}	
 

double tot_func(double *x, double *p)
{
  if (x[0]<m_min || x[0]>m_max) return 0.;
  return sgn_func(x,p)+ bkg_func(x,&p[3]);
}



void fit_mc()
{ 
  //TFile * f = new TFile("check_here_for_the_plots_in_charmWG.root");

  TFile * f = new TFile("histo_ksks_proc11_bucket_9to12_unscaled.root");  
  //TFile * f = new TFile("histo_ksks_proc11_unscaled.root");
  //TFile * f = new TFile("histo_ksks_proc11_scaled.root");
  //TFile * f = new TFile("histo_ksks_bucket_9to12_unscaled.root");
  //TFile * f = new TFile("histo_ksks_bucket_9to12_scaled.root");
  //TFile * f = new TFile("histo_ksks_proc11_bucket_9to12_scaled.root");

  TH1F *h1 = new TH1F("h1","histogram",100, 0.139, 0.159);
   
  // h1 = (TH1F*)f-> Get("DeltaM_fde");
  h1 = (TH1F*)f-> Get("DeltaM4");
  
  gStyle->SetOptFit(1111);
  gStyle->SetHistMinimumZero(kTRUE);
  
  // Fit
  TF1 *tot = new TF1("tot",tot_func,m_min,m_max,7);
  tot->SetNpx(n_bins*3);
    
  //tot->SetParameters(2e3,0.145,0.0006,1e2,9,20,0.1399);//ft >0.001 to 0.003
  tot->SetParameters(2e3,0.145,0.0006,1e2,9,20,0.1399);//fte > 4
 
  tot->SetParNames("N_{sig}","#mu","#sigma","N_{bkg}","#alpha","#beta","m0");
  
  tot->SetParLimits(1,m_min,m_max);//mu1
  tot->SetParLimits(2,0.,10.);//sigma1
  
  tot->FixParameter(6,0.1399); //m0
  tot->FixParameter(5,0.); //beta
  
  //h1->Fit(tot,"R");
  h1->Fit(tot,"L");
  h1->SetMarkerStyle(20);
  h1->Draw("E");
  
  TF1 *sgn = new TF1("sgn",sgn_func,m_min,m_max,3);
  sgn->SetLineColor(kGreen);
  sgn->SetLineStyle(2);
  sgn->SetParameters(tot->GetParameter(0),tot->GetParameter(1),tot->GetParameter(2));
  //sgn->SetParameters(tot->GetParameters());
  sgn->Draw("same");
  
  TF1 *bkg = new TF1("bkg",bkg_func,0.139,0.159,4);
  bkg->SetLineColor(kBlue);
  bkg->SetLineStyle(2);
  bkg->SetParameters(tot->GetParameter(3),tot->GetParameter(4),tot->GetParameter(5),tot->GetParameter(6));
  bkg->Draw("same");
    
  double w=h1->GetBinWidth(1);
  
  double integral0 = sgn->Integral(0.1445,0.1465)/w;
  double integral1 = bkg->Integral(0.1445,0.1465)/w;
  
  cout<<"Integral of sgn in [0.1445,0.1465]= "<< integral0 <<endl;
  cout<<"Integral of bkg in [0.1445,0.1465]= "<< integral1 <<endl;


}