const unsigned int n_bins(100);
const double m_min(1.8), m_max(1.95);
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 f1 = p[1];
    double mu1 = p[2];
    double sigma1 = p[3];
    double mu2 = p[4];
    double sigma2 = p[5];
    return Nsig*(f1*gaus_func(m,mu1,sigma1)+(1.-f1)*gaus_func(m,mu2,sigma2))*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 lambda = p[1];
   if (lambda==0.) return 1./(m_max-m_min);
   double norm = (TMath::Exp(-lambda*m_min)-TMath::Exp(-lambda*m_max))/lambda;
   if (norm==0.) return 0.;
   return Nbkg*TMath::Exp(-lambda*m)/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[6]);
}





  void binned_D01_Kpipi0_RS_M_D0()
{

   TFile * f = new TFile("D0_M_all_set_kpipi0_total_proc11_4loop.root");

   TFile * f1 = new TFile("D0_M_all_set_kpipi0_total_proc11_4loop_after_fitting.root","recreate");


   int histocounter = 0;
   //Creat an iterator
   TIter nextkey (f->GetListOfKeys());
   TKey *key;

   //Loop over the keys in the file
   while ((key = (TKey*) nextkey())) {

      //Read object from first source file
      TObject *obj = key ->ReadObj();

      TH1* h1 = (TH1*) (obj);
      h1-> Print();

      gStyle->SetOptFit(1111);
      gStyle->SetHistMinimumZero(kTRUE);

      // Fit
      TF1 *tot = new TF1("tot",tot_func,m_min,m_max,8);
      tot->SetNpx(n_bins*3);

      tot->SetParameters(5e5,0.7,1.85,0.01,1.87,0.007,2e1,1.);  // jan recom with pCMS 2 & without Dz_pt

      tot->SetParNames("N_{sig}","f_{1}","#mu_{1}","#sigma_{1}","#mu_{2}","#sigma_{2}","N_{bkg}","#lambda");
      tot->SetParLimits(1,0.,1.);//f1
      tot->SetParLimits(2,m_min,m_max);//mu1
      tot->SetParLimits(3,0.,10.);//sigma1
      tot->SetParLimits(4,m_min,m_max);//mu2
      tot->SetParLimits(5,0.,10.);//sigma2

      h1->Fit(tot,"R");
      h1->SetMarkerStyle(20);
      h1->Draw("E");

      TF1 *sgn = new TF1("sgn","sgn_func",m_min,m_max,6);
      sgn->SetLineColor(kGreen);
      sgn->SetLineStyle(2);
      sgn->SetParameters(tot->GetParameter(0),tot->GetParameter(1),tot->GetParameter(2),tot->GetParameter(3),tot->GetParameter(4),tot->GetParameter(5));
      sgn->Draw("same");

      // Draw bkg projection
      TF1 *bkg = new TF1("bkg","bkg_func",m_min,m_max,2);
      bkg->SetLineColor(kBlue);
      bkg->SetLineStyle(2);
      bkg->SetParameters(tot->GetParameter(6),tot->GetParameter(7));
      bkg->Draw("same");

      gPad->Update();
      histocounter++;
      gPad->SaveAs(Form("hist%d.pdf",histocounter));

      // efficiency and purity

      double w= h1->GetBinWidth(1);
      double integral0 = sgn->Integral(1.8,1.95)/w;
      cout<<"Integral of sgn in [1.8, 1.95]= "<< integral0 <<endl;

      double Entries = h1->GetEntries();
      cout<< "Total Entries =" << Entries<< endl;

      double purity = (integral0/Entries)* 100;
      cout<<"purity in % ="<< purity << endl;

      cout<<"Efficiency w.r.t ICHEP results ="<< integral0/30680 << endl;
      cout<<"purity w.r.t ICHEP results ="<< purity/88.36 << endl;
      cout<< "==========================================================================================" << endl;

      h1->Write();
   }

   f1->Close();

}