const unsigned int n_bins(100);
const double m_min(0.1396), 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 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 mu1 = p[3];
  double sigma1 = p[4];
  double gamma = p[5];
  double delta = p[6];
  double mu2 = p[7];
  double sigma2 = p[8];
  double mu3 = p[9];
  double sigma3 = p[10];
  return Nsig*(f1*johnson(m,mu1,sigma1,gamma,delta) + (1.-f1)*f2*gaus_func(m,mu2,sigma2)+ (1.-f1)*(1-f2)*gaus_func(m,mu3,sigma3))*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[11]);
}

void kk_fit_full_set()
{          
  TFile * f = new TFile("deltam_dist_all_pid_comb.root");
  TFile * f1 = new TFile("deltam_full_set_kk_pid_comb.root","RECREATE");
  
  TCanvas canvas("canvas");
  canvas.Print("deltam_full_set.pdf[");
  
  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);
    
    TF1 *tot = new TF1("tot",tot_func,m_min,m_max,15);
    tot->SetNpx(n_bins*3);
    
    //tot->SetParameters(7e6, 0.4, 0.2, 0.145, 0.002, 1., 5., 0.146, 0.002, 0.146, 0.001);//data (don't use these param)
    tot->SetParameters(7e6, 0.5, 0.2, 0.145, 0.002, 1., 5., 0.146, 0.002, 0.146, 0.001);//mc
    //tot->SetParameters(6e6, 0.5, 0.2, 0.145, 0.002, 1., 5., 0.1445, 0.002, 0.144, 0.001);//mc (to fix PID(k, pi) 0.5)
    tot->SetParameter(11,5e4);//mc
    //tot->SetParameter(11,7e5);//data
    tot->SetParameter(12,8);
    tot->SetParameter(13,12);
    tot->SetParameter(14,0.1396);  
   
    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, "#alpha");
    tot->SetParName(13, "#beta");
    tot->SetParName(14, "m0");
    
    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
    
    tot->FixParameter(13,0); //beta
    tot->FixParameter(14,0.1396); //m0
    
    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->SetParameters(tot->GetParameters());
    sgn->Draw("same");
    
    
    TF1 *bkg = new TF1("bkg",bkg_func,0.1396,0.159,4);
    bkg->SetLineColor(kBlue);
    bkg->SetLineStyle(2);
    bkg->SetParameters(tot->GetParameter(11),tot->GetParameter(12),tot->GetParameter(13),tot->GetParameter(14));
    bkg->Draw("same");

        
    canvas.Clear();
    h1->Draw();
    canvas.Print("deltam_full_set.pdf");
     //calculations
    
    double w= h1->GetBinWidth(1);
    double integral0 = sgn->Integral(0.139,0.159)/w;
    double integral1 = bkg->Integral(0.139,0.159)/w;
    double integral2 = sgn->Integral(0.1445,0.1465)/w;
    double integral3 = bkg->Integral(0.1445,0.1465)/w;

    cout<<"Integral of sgn in full region= "<< integral0 <<endl;
    cout<<"Integral of bkg in full region= "<< integral1 <<endl;
    cout<<"Integral of sgn in [0.1445, 0.1465]= "<< integral2 <<endl;
    cout<<"Integral of bkg in [0.1445, 0.1465]= "<< integral3 <<endl;
    //cout<<"Signal Yield = "<< N_{sig} <<endl;
    
       
    
    h1->Write();
  }
  canvas.Print("deltam_full_set.pdf]");   // Finished; close the multipage file
  f1->Close();
}