const unsigned int n_bins(100);
const double m_min(0.45), m_max(0.55);
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 johnson(double m, double mu, double sigma, double gamma, double delta)
{
  if(m<m_min || m>m_max) return 0.;
  if(sigma==0. || delta==0.) return 0.;
  double arg_max = gamma+delta*TMath::ASinH((m_max-mu)/sigma);
  double arg_min = gamma+delta*TMath::ASinH((m_min-mu)/sigma);
  double norm    = TMath::Sqrt(TMath::Pi()/2.)*sigma/delta*(TMath::Erf(arg_max/TMath::Sqrt2())-TMath::Erf(arg_min/TMath::Sqrt2()));
  if(norm<=0.) return 0.;
  double z     = (m-mu)/sigma;
  double arg   = gamma+delta*TMath::ASinH(z);
  double value = TMath::Exp(-arg*arg/2.)/TMath::Sqrt(1.+z*z);
  return (value>0.) ? value/norm : 0.;
}

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 f2 = p[2];
    double mu_J = p[3];
    double sigma_J = p[4];
    double gamma = p[5];
    double delta = p[6];
    double mu1 = p[7];
    double sigma1 = p[8];
    double mu2 = p[9];
    double sigma2 = p[10];
    
    return Nsig*(f1*johnson(m,mu_J,sigma_J,gamma,delta) + (1.-f1)*f2*gaus_func(m,mu1,sigma1)+ (1.-f1)*(1-f2)*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[11]);//till p10, it is signal
}

void binned_fit_KS_mass()
{          
  TFile * f = new TFile("KS_pipi_histo_b17_05_02_05_200bin.root");   
  TH1 *h1 = new TH1F("h1","histogram",100, 0.45, 0.55);
  h1 = (TH1*)f-> Get("Ks_7");
  //h1 = (TH1F*)f-> Get("Ks_6");
  
  gStyle->SetOptFit(1111);
  gStyle->SetHistMinimumZero(kTRUE);

  TF1 *tot = new TF1("tot",tot_func,m_min,m_max,13); //sig:11 and bkg:2
  tot->SetNpx(n_bins*3);

  //100 bins
  //tot->SetParameters(2e6, 0.99, 0.43, 0.49, 0.002, 1., 5., 0.49, 0.002, 0.488, 0.001);
  //tot->SetParameter(11,8e5);
  //tot->SetParameter(12,-9);

  //F0R  BASELINE(200 bins)
  tot->SetParameters(2e6, 0.9, 0.56, 0.49, 0.002, 1., 5., 0.49, 0.002, 0.488, 0.001);
  tot->SetParameter(11,13e5);
  tot->SetParameter(12,-1);

  
  //F0R  BASELINE(500 bins)
  //tot->SetParameters(2e6, 0.99, 0.41, 0.49, 0.002, 1., 5., 0.49, 0.002, 0.488, 0.001);
  //tot->SetParameter(11,15e5);
  //tot->SetParameter(12,-2); 
    
  tot->SetParNames("N_{sig}", "f_{J}", "f_{G1}", "#mu_{J}", "#sigma_{J}", "#gamma_{J}", "#delta_{J}", "#mu_{G1}", "#sigma_{G1}", "#mu_{G2}", "#sigma_{G2}");


  //F0R  ALTERNATIVE
  //tot->SetParameters(2e6, 0.99, 0.41, 0.49, 0.002, 1., 5., 0.49, 0.002, 0.488, 0.001);
  //tot->SetParameter(11,15e5);
  //tot->SetParameter(12,-10);   
  //tot->SetParNames("N_{sig}", "f_{J}", "f_{G1}", "#mu_{J}", "#sigma_{J}", "#gamma_{J}", "#delta_{J}", "#mu_{G1}", "#sigma_{G1}", "#mu_{G2}", "#sigma_{G2}");
  tot->SetParName(11, "N_{bkg}");
  tot->SetParName(12, "#lambda");
  tot->SetParLimits(1,0.,1.);//f1
  tot->SetParLimits(2,0.,1.);//f1
  tot->SetParLimits(3,m_min,m_max);//mu1
  tot->SetParLimits(4,0.,10.);//sigma1
  tot->SetParLimits(7,m_min,m_max);//mu2
  tot->SetParLimits(8,0.,10.);//sigma2
  tot->SetParLimits(9,m_min,m_max);//mu2
  tot->SetParLimits(10,0.,10.);//sigma2
  
  //Pull Plot
  TH1 *hpull = new TH1F("hpull","",200, 0.45, 0.55);
  Int_t bins = h1->GetNbinsX();
  cout << "no. of bins = " << bins << endl;

  cout<< h1->GetBinContent(10)<<endl;
  cout<< tot->Eval(h1->GetBinCenter(10))<<endl;
  cout<<h1->GetBinContent(10) - tot->Eval(h1->GetBinCenter(10)) <<endl;
  
  for (int i =1;i<=bins;i++)
    {
      Double_t pull = (h1->GetBinContent(i) - tot->Eval(h1->GetBinCenter(i)));
      cout << "pulls = " << pull << endl;
      hpull->Fill(pull);
    }
  TCanvas *canvas = new TCanvas("canvas","canvas", 600, 600);

  Double_t xlow, ylow, xup, yup;
  canvas->GetPad(0)->GetPadPar(xlow,ylow,xup,yup);
  canvas->Divide(1,2);

  TVirtualPad *upPad = canvas->GetPad(1);
  upPad->SetPad(xlow,ylow+0.25*(yup-ylow),xup,yup);
  
  TVirtualPad *dwPad = canvas->GetPad(2);
  dwPad->SetPad(xlow,ylow,xup,ylow+0.25*(yup-ylow));

  //canvas->Update();
  canvas->cd(1);
  
  //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,11);
  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),tot->GetParameter(6),tot->GetParameter(7),tot->GetParameter(8),tot->GetParameter(9), tot->GetParameter(10));
  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(11),tot->GetParameter(12));
  bkg->Draw("same");
  
  canvas->cd(2);
  hpull->Draw();
  
  canvas->Update();
   
}