// Normalizing the Gaussian fit - subtracting off the exponential portion.

//#include <TAxis.h>

Int_t npoints;
Double_t norm = 1;

const double EXPO_SLOPE = 1.9;
const double EXPO_MAXSLOPE = 2.0;
const double EXPO_MINSLOPE = 1.8;
const double MASS_JPSI = 3.097;
const double Max_bins = 320;

Double_t xmin, xmax;
Int_t ibins;

//Crystal ball function added June 28
double CBFunction(double *x, double *p)// crystal ball function
{
  double norm = p[0];
  double alpha = p[1];  // tail start 1.0
  double n = p[2];      // tail shape  0.7
  double mu = p[3];     //J/Psi mass 3.0967
  double sigma_JPsi = p[4];  //J/Psi width 0.075

  double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
  double B = n/fabs(alpha) - fabs(alpha);
  double k = (x[0]-mu)/sigma_JPsi ;
  double val;
  if( k > -alpha )
  {
    val = norm*TMath::Exp(-0.5*pow(k,2));//val =  Norm * exp( -pow(x-x_mean,2) / (2.0*pow(sigma,2)));
  // cout << " x " << " x_mean " << x_mean << " sigma " << sigma << " Norm " << Norm << " val " << val << endl;
  }
  else
  {
    val = norm*A*pow(B-k,-n);//val = Norm * A * pow(B - (x-x_mean)/sigma, -n);
  }

}

double CBFunction_bkg(double *x, double *p) // // added an exponential background to Crystal ball function July 1,2019
{
  double norm = p[0];
  double alpha = p[1];  // tail start 1.0
  double n = p[2];      // tail shape  0.7
  double mu = p[3];     //J/Psi mass 3.0967
  double sigma = p[4];  //J/Psi width 0.075

  double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
  double B = n/fabs(alpha) - fabs(alpha);
  double k = (x[0]-mu)/sigma;
  double val;
  if( k > -alpha )
  {
    val = norm*TMath::Exp(-0.5*pow(k,2));//val =  Norm * exp( -pow(x-x_mean,2) / (2.0*pow(sigma,2)));
  // cout << " x " << " x_mean " << x_mean << " sigma " << sigma << " Norm " << Norm << " val " << val << endl;
  }
  else
  {
    val = norm*A*pow(B-k,-n);//val = Norm * A * pow(B - (x-x_mean)/sigma, -n);
  }

  double bkgnorm1 = p[5];//Norm of the exponential
  double bkgslope = p[6]; //slope  of the background

  double bkg = bkgnorm1 * exp(-bkgslope * x[0]);

  val = val + bkg;
  
  if( TMath::IsNaN(val) )
  {
    val = 0.0;
  }

  return val;
}


//======================================================================
//Reading in the file, cloning the histogram.
//======================================================================

void InitialCBnormfit(const char fname = "hnorm.root")
{
  //gStyle->SetOptFit(111111);
  gStyle->SetOptStat(kFALSE);

  //Get file and get  histogram
  TFile *file = TFile::Open("hnorm.root");
  TH1F *hmass = (TH1F*)file->Get("hmass_copy");
  TH1F *cbhmass = (TH1F*)file->Get("CBh_mass");
  TH1F *hsub = (TH1F*)file ->Get("hmass_copy")->Clone("hsub");


  TH1F *hsubmin = (TH1F*)file ->Get("hmass_copy")->Clone("hsubmin");
  TH1F *hsubmax = (TH1F*)file ->Get("hmass_copy")->Clone("hsubmax");

  //TH1F *hsub = (TH1F*)file ->Get("hmass_copy")->Clone("mass_sub");

  hsub->Reset();
  hsubmin->Reset();
  hsubmax->Reset();
  hsubmax->Sumw2();
  hsubmin->Sumw2();
  cbhmass->Sumw2();


  TH1F *hdndm = (TH1F*)file ->Get("hmass_copy")->Clone("hdndm");
  //TH1F *CBh_dndm = (TH1F*)file ->Get("hmass_copy")->Clone("CBh_dndm");//Crystal ball histogram to fit crystal ball - added June 28
  TH1F *CBh_dndm = (TH1F*)file ->Get("CBh_mass")->Clone("CBh_dndm");



  hdndm->Sumw2();
  CBh_dndm->Sumw2();
  int bin = 1;


  //====================================================================
  //Canvas dndm
  //===================================================================
  TCanvas *ac = new TCanvas("cdndm","dN/dm",550,425);
  TCanvas *bc = new TCanvas("cbdndm","cbdN/dm",550,425);


  //===================================================================
  // Calculating the counts, normalizing and plotting dn/dm
  //===================================================================
  Double_t binwid = hdndm->GetBinWidth(1);
  Double_t scale = 1.0/(hdndm-> GetBinWidth(1));//Integral(bmin,bmax); 
  hdndm->Scale(scale);
  hdndm->Draw("ehist");
  Double_t integ = hdndm->Integral("width");
  cout << "dN/dm integ " << integ << " = " << hmass_copy->Integral() << endl;



//====================================================================
  //calculating the J/Psi count from crystal ball fit
 //====================================================================
  Double_t binwidth = CBh_dndm->GetBinWidth(1);
  Double_t scale = 1.0/(CBh_dndm-> GetBinWidth(1));//Integral(bmin,bmax); 
  hdndm->Scale(scale);
  hdndm->Draw("ehist");
  Double_t CBinteg = CBh_dndm->Integral("width");
  cout << "dN/dm CB-integ " << CBinteg << " = " << CBh_mass->Integral() << endl;



  //====================================================================
  // Fitting normalized Gaussian + exponential to the function.
  // Fixing the J/Psi mass and the slope of the exponential. 
  //====================================================================

  TF1* fit = new TF1 ("fit", "[0]*exp(-0.5*pow((x-[1])/[2],2))/([2]*sqrt(2 *TMath::Pi())) + [3]*exp(-[4]*x)", 2.4, 4.0);
  fit->SetParNames("amp_gauss", "Mean_Gaus", "sigma_Gaus", "amp_exp", "exp_arg");
  fit->SetParameters(4, 3.1, 0.1, 1000, 1.7);
  fit->FixParameter(1,MASS_JPSI);
  fit->FixParameter(4,EXPO_SLOPE);
  //hdndm->Fit("fit", "WL", "", 2.41, 4.0);
  hdndm->Fit(fit, "LLI", "", 2.41, 4.0);//not printing out on canvas, why?
  hdndm->GetXaxis()->SetRangeUser(2.4,4.0);
  fit->SetLineColor(kRed);
  hdndm->Draw("ehist");
  fit->Draw("same");



 //===========================================================================================  
 //fitting Guassian only to find counts: July 3, 2019
//============================================================================================

  TF1 *normfit = new TF1("normfit", "[0]*exp(-0.5*pow((x-[1])/[2],2))/([2]*sqrt(2 *TMath::Pi()))",2.4, 4.0);
  normfit->SetParNames("amp_gauss", "Mean_Gaus", "sigma_Gaus");
  normfit->SetParameters(1.631,3.097,0.05718);
  normfit->SetParameter(0,fit->GetParameter(0));

  hdndm->Fit(normfit, "RL", "", 2.4, 3.2);


  //-------------------------------------------------------------------------
  // Fitting normalized crystal ball + exponential.
  //  Fixing exponential slope and mass of J/Psi  added July 2019
  // ------------------------------------------------------------------------

 
  TF1* cbfit = new TF1("cbfit",CBFunction_bkg,2.4,5.,7);
  cbfit->SetLineColor(kBlue);
  cbfit->SetParNames("norm","alpha","n","mu","sigma","bkgnorm","bkgslope");
  cbfit->SetParameters(16,1.0,1.5.,3.097,0.07,100.,2.);
  //cbfit->FixParameter(3,MASS_JPSI);
 // cbfit->FixParameter(6,EXPO_SLOPE);
  CBh_dndm->GetXaxis()->SetTitle("mass(GeV/c^2)");
  CBh_dndm->GetYaxis()->SetTitle("Counts");
  CBh_dndm->Fit("cbfit","RL","",2.4,4.0);
  CBh_dndm->Draw("ehist");
  CBh_dndm->GetXaxis()->SetRangeUser(2.4,4.0);
  cbfit->SetLineColor(kMagenta);
  CBh_dndm->Draw("ehist");
  cbfit->Draw("same");


  TF1 *cnormfit = new TF1("cnormfit", "[0]*exp(-0.5*pow((x-[1])/[2],2))/([2]*sqrt(2 *TMath::Pi()))",2.4, 4.0);
  cnormfit->SetLineColor(kRed);
  cnormfit->SetParNames("amp_gauss", "Mean_Gaus", "sigma_Gaus");
  cnormfit->SetParameters(1.631,3.097,0.05718);
  cnormfit->SetParameter(0,fit->GetParameter(0)); 
  CBh_dndm->Fit("cnormfit","RL","",2.4,3.2);
  cnormfit->Draw("same");




  //=====================================================================
  //getting the normalization parameter, error, mean and the amplitude
  // for Gaussian + exponential and the crystal ball + exponential 
  //======================================================================
  Double_t njpsi = fit->GetParameter(0);
  Double_t njpsi_err = fit->GetParError(0);
  Double_t sigma = fit->GetParameter(1);
  Double_t sigma2 = fit->GetParameter(2);

  Double_t bkg = fit->GetParameter(3);
  Double_t bkg_err = fit->GetParError(3);
  // bkg *= hsub->GetBinWidth(1); bkg_err *= hsub->GetBinWidth(1); // same scaling
  cout << "Fit N_J/Psi  " << njpsi << "\t" << njpsi_err << endl;


 // CBh_dndm->Draw("e");
  Double_t njpsi_cb = cbfit->GetParameter(0);
  Double_t njpsi_errcb = cbfit->GetParError(0);
  Double_t sigmacb = cbfit->GetParameter(3);
  Double_t sigmacb2 = cbfit->GetParameter(4);

  Double_t cbbkg = cbfit->GetParameter(5);
  Double_t cbbkg_err = cbfit->GetParError(5);
  // bkg *= hsub->GetBinWidth(1); bkg_err *= hsub->GetBinWidth(1); // same scaling
  cout << "Fit CBN_J/Psi  " << njpsi_cb << "\t" << njpsi_errcb << endl;

    

  //====================================================================
  //subtracting exponential from the normalized Gaussian function
  //====================================================================


  TCanvas *cc = new TCanvas("ccounts","counts",800,800);
  cc->Divide(2,2);

  cc->cd(1);
  TF1* fit2 = new TF1 ("fit2", "[0]*exp(-0.5*pow((x-[1])/[2],2))/([2]*sqrt(2 *TMath::Pi())) + [3]*exp(-[4]*x)", 0., 1000.);
  fit2->SetParNames("amp_gauss", "Mean_Gaus", "sigma_Gaus", "amp_exp", "exp_arg");
  fit2->SetParameters(4, 3.1, 0.1, 1000,1.7);
  fit2->FixParameter(1,MASS_JPSI);
  fit2->FixParameter(4,EXPO_SLOPE);
  hmass->Fit("fit2", "L", "", 2.4, 4.);

  TF1 * fitexpo = new TF1("fitexpo","[0]*exp(-[1]*x)", 0., 1000.);//exponential function
  fitexpo->SetParameters(fit2->GetParameter(3), EXPO_SLOPE);
  fitexpo->SetLineColor(kBlue);//Draw using different color

  TF1 * CBfitexpo = new TF1("fitexpo","[0]*exp(-[1]*x)", 0., 1000.);//exponential function
  CBfitexpo->SetParameters(fit2->GetParameter(3), EXPO_SLOPE);
  CBfitexpo->SetLineColor(kMagenta);//Draw using different color

 // cout << "*** fit to counts ***" << endl;
 // cout << "expo bkg " << fit2->GetParameter(3) << "\t" << endl;

  TF1 * fmaxexpo = new TF1("fmaxexpo","[0]*exp(-[1]*x)", 0., 1000.);//exponential function
  fmaxexpo->SetLineColor(kGreen);//Draw using different color
  TF1 * fminexpo = new TF1("fminexpo","[0]*exp(-[1]*x)", 0., 1000.);//exponential function
  fminexpo->SetParameters(fit2->GetParameter(3), EXPO_MINSLOPE);
  fminexpo->SetLineColor(kYellow);//Draw using different color



  //======================================================================================
    //subtracting exponential from the Crystal ball function added July 2, 2019: causing crash in root,
  //======================================================================================


  TCanvas *fc = new TCanvas("ccounts","counts",800,800);
  fc->Divide(2,2);

  fc->cd(1);
  TF1* cbfit2 = new TF1 ("cbfit2", CBFunction_bkg, 2.4,5.,7);
  cbfit->SetParNames("norm","alpha","n","mu","sigma","bkgnorm","bkgslope");
  cbfit->SetParameters(16,1.0,1.5.,3.097,0.07,100.,2.);
  cbfit2->FixParameter(3,MASS_JPSI);
  cbfit2->FixParameter(6,EXPO_SLOPE);
  CBh_mass->Fit("cbfit2", "L", "", 2.4, 4.);

  TF1 * fitexpo_CB = new TF1("fitexpo_CB","[0]*exp(-[1]*x)", 0., 1000.);//exponential function
  fitexpo_CB->SetParameters(cbfit2->GetParameter(5), EXPO_SLOPE);
  fitexpo_CB->SetLineColor(kBlue);//Draw using different color

   TF1 * Cfitexpo = new TF1("fitexpo_CB","[0]*exp(-[1]*x)", 0., 1000.);//exponential function
  Cfitexpo->SetParameters(cbfit2->GetParameter(5), EXPO_SLOPE);
  Cfitexpo->SetLineColor(kMagenta);//Draw using different color

  //cout << "*** fit to counts ***" << endl;
 // cout << "expo bkg " << cbfit2->GetParameter(5) << "\t" << endl;

  TF1 * cbfmaxexpo = new TF1("cbfmaxexpo","[0]*exp(-[1]*x)", 0., 1000.);//exponential function
  cbfmaxexpo->SetLineColor(kGreen);//Draw using different color
  TF1 * cbfminexpo = new TF1("cbfminexpo","[0]*exp(-[1]*x)", 0., 1000.);//exponential function
  cbfminexpo->SetParameters(cbfit2->GetParameter(5), EXPO_MINSLOPE);
  cbfminexpo->SetLineColor(kYellow);//Draw using different color

  CBh_mass->Draw();
  CBh_mass->GetXaxis()->SetRangeUser(2.5,4.0);
  fitexpo_CB->DrawCopy("same");




  //======================================================================================
  // Drawing the backgrounds on the hmass
  // =====================================================================================
  hmass->Draw();
  hmass->GetXaxis()->SetRangeUser(2.5,4.0);
  fitexpo->DrawCopy("same");
  

  //=======================================================================================
  // fixing of the min and max background of the normal Gaussian
  // Fixing background for the Crystal ball added July 2019
  //=======================================================================================
 
  fmaxexpo->SetParameters(fit2->GetParameter(3), EXPO_MAXSLOPE);
  fmaxexpo->DrawCopy("same");

  fminexpo->SetParameters(fit2->GetParameter(3), EXPO_MINSLOPE);
  fminexpo->DrawCopy("same");

  cc->cd(2);
  int nbins = hmass->GetNbinsX();
  for (int ibin=1; ibin<=nbins; ibin++)
  {
    double x = hmass->GetXaxis()->GetBinCenter(ibin);
    double fval = fitexpo->Eval(x);
    hsub->SetBinContent(ibin,hmass->GetBinContent(ibin)-fval);
    /*
       cout << ibin << "\t hmass = " << hmass->GetBinContent(ibin) << "\tfval = " << fval  <<"\t hsub = "
       << hsub->GetBinContent(ibin) << endl;
       */
  }  
  
//not working yet.
  cbmaxexpo->SetParameters(fit2->GetParameter(5), EXPO_MAXSLOPE);
  cbmaxexpo->DrawCopy("same");

  cbminexpo->SetParameters(fit2->GetParameter(5), EXPO_MINSLOPE);
  cbminexpo->DrawCopy("same");

  fc->cd(2);
  int nbins = CBh_mass->GetNbinsX();
  for (int ibin=1; ibin<=nbins; ibin++)
  {
    double x = CBh_mass->GetXaxis()->GetBinCenter(ibin);
    double cbval = fitexpo_CB->Cval(x);
    chsub->SetBinContent(ibin,CBh_mass->GetBinContent(ibin)-cbal);
    
  }

   

  //=============================================================================
  // plot regular background subtracted mass, plotting the e+e- continuum 
  // that is +- 1 sigma of the background
  // ============================================================================
  hsub->Draw("e");
  hsub->GetXaxis()->SetRangeUser(2.4,4.0);
  TLine *line = new TLine(2.4,0,4.0,0);
  line->SetLineColor(kRed);
  line->Draw();

  // +1 sigma
  cc->cd(3);
  int nbins = hmass->GetNbinsX();
  for (int ibin=1; ibin<=nbins; ibin++)
  {
    double x = hmass->GetXaxis()->GetBinCenter(ibin);
    double fval = fmaxexpo->Eval(x);
    hsubmax->SetBinContent(ibin,hmass->GetBinContent(ibin)-fval);
    /*
       cout << ibin << "\t hmass = " << hmass->GetBinContent(ibin) << "\tfval = " << fval  <<"\t hsub = "
       << hsub->GetBinContent(ibin) << endl;
       */
  }   

  hsubmax->Draw("E1");
  hsubmax->GetXaxis()->SetRangeUser(2.4,4.0);

  // -1 sigma
  cc->cd(4);
  int nbins = hmass->GetNbinsX();
  for (int ibin=1; ibin<=nbins; ibin++)
  {
    double x = hmass->GetXaxis()->GetBinCenter(ibin);
    double fval = fminexpo->Eval(x);
    hsubmin->SetBinContent(ibin,hmass->GetBinContent(ibin)-fval);
  
  }   

  hsubmin->Draw("E1");
  hsubmin->GetXaxis()->SetRangeUser(2.4,4.0);


  //=========================================================================
  // Counting the number of  J/Psi in the sample
  // ======================================================================== 

  float start_int = 2.81;
  float stop_int = 3.39;


  int fbin = hmass->FindBin(start_int);
  int lbin = hmass->FindBin(stop_int);
  double sum = 0;
  double sum_integ = 0;
  double sumerr = 0;


  for(int k=fbin; k<=lbin; k++)
  {
    float bincenter = hmass->GetBinCenter(k);
    sum += (hmass->GetBinContent(k) - fitexpo->Eval(bincenter));
    sum_integ +=(fit->Eval(bincenter) - fitexpo->Eval(bincenter));
    sumerr += pow(hmass->GetBinError(k),2);
  }

  cout<<"Jpsi direct count in mass range 2.8 < m_ee < 3.4 =" << sum <<" +- "<< sqrt(sumerr)<<endl;

  //cout <<"Jpsi raw data counts = " <<hmass->Integral(2.8,3.4)<<endl;


  //==================================================================
  //analyzing the error in the e+e- continuum
  //==================================================================


  double sum_min = 0;
  double sum_integ = 0;
  double sumer_min = 0;


  int bbins = hsubmin-> FindBin(start_int);
  int ubins = hsubmin->FindBin(stop_int);
  for(int l=bbins; l<=ubins; l++)
  {
    float bincenter = hsubmin->GetBinCenter(l);
    sum_min += hsubmin->GetBinContent(l);
    //sum_integ +=(fit->Eval(bincenter) - fminexpo->Eval(bincenter));
    sumer_min += pow(hsubmin->GetBinError(l),2);

    // cout<<"min: "<<l<<" "<<sum_min <<endl;
    // cout << "min_sub = " << hsubmin->GetBinContent(l) << "  fminex = " << fminexpo->Eval(bincenter)<<endl;
  }
  cout<<"Background error min = " << sum_min <<" +- "<< sqrt(sumer_min)<<endl;


  double sum_max = 0;
  double sum_integ = 0;
  double sumer_max = 0;


  int b_bins = hsubmax-> FindBin(start_int);
  int u_bins = hsubmax->FindBin(stop_int);


  for(int j =b_bins; j<=u_bins; j++)
  {
    float bincenter = hsubmax->GetBinCenter(j);
    sum_max += hsubmax->GetBinContent(j);
    //sum_integ +=(fitexpo->Eval(bincenter) - fmaxexpo->Eval(bincenter));
    sumer_max += pow(hsubmax->GetBinError(j),2);
  }
  cout<<"Background error max = " << sum_max <<" +- "<< sqrt(sumer_max)<<endl;

  //cout<< "Difference = "<<sum_max -sum_min<<endl; 

  //cout <<"Bkg min error  = " <<hsubmin->Integral(2.6,3.8)<<endl;

  //cout <<"Bkg max error  = " <<hsubmax->Integral(2.6,3.8)<<endl;
}