//
//  mass plots and background.cpp
.
//

#include "mass plots and continuum background.hpp"

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

//#include <TAxis.h>

Int_t npoints,bin;
Double_t norm = 1;

void Initalnormfit(const char fname = "hnorm.root")
{
    //Get file and get  histogram
    TFile *file = TFile::Open("hnorm.root");
    TH1F *hmass = (TH1F*)file ->Get("hmass_copy");
    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");
    
    hsub->Reset();
    hsubmin->Reset();
    hsubmax->Reset();
    
    TH1F *hdndm = (TH1F*)file ->Get("hmass_copy")->Clone("hdndm");
    hdndm->Sumw2();
    
    //TH1F * f = new TH1F
    
    //hnorm-> GetBinContent(bin);
    int bin = 1;
    //Double_t x,xmin,xmax;
    //int bin = hnorm->GetXaxis()->FindBin(x);
    //int bmin = GetXaxis->FindBin(xmin);
    //int bmax = GetXaxis->FindBin(xmax);
    
    
    Double_t binwid = hdndm->GetBinWidth(1);
    Double_t scale = 1.0/(hdndm-> GetBinWidth(1));//Integral
    hdndm->Scale(scale);
    hdndm->Draw("e");
    Double_t integ = hdndm->Integral("width");
    cout << "dN/dm integ " << integ << " = " << hmass->Integral() << endl;
    
    TF1* fit = new TF1 ("fit", "[0]*exp(-0.5*pow((x-[1])/[2],2))/([2]*sqrt(2 *TMath::Pi())) + [3]*exp(-[4]*x)", 0., 1000.);

    fit->SetParNames("amp_gauss", "Mean_Gaus", "sigma_Gaus", "amp_exp", "exp_arg");
    fit->SetParameters(60, 3.1, 0.2, 1000, 4.0);
    //fit->FixParameter(3.097, mean_Gaus);?
    hdndm->Fit("fit", "WL", "", 2.6, 4.);
    
    //getting the normalization parameter, error, mean and the amplitude.
    Double_t jpsi_n = fit->GetParameter(0);
    Double_t jpsi_errn = 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
    
    //Double_t bkg = fit->GetParameter(3)*hmass->GetBinWidth(1);
    //Double_t bkg_err = fit->GetParError(3);
    cout << "N_J/Psi  " << jpsi_n << "\t" << jpsi_errn << endl;
    
    cout << "bkg = " << bkg << " \t" << "bkg_err = " << bkg_err << endl;
    
    //Canvas  for dndm
    TCanvas *ac = new TCanvas("cdndm","dN/dm",550,425);
    //TCanvas *c1 = new TCanvas("c1","multipads",900,700);
    hdndm->Draw();
    hdndm->GetXaxis()->SetRangeUser(2.6,4.0);
    fit->SetLineColor(kRed);
    fit->Draw("same");
    
    double sigma;
    
    
    //==============================
    
    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(20, 3.1, 0.2, 1000, 4.0);
    //fit2->FixParameter(3.097, mean_Gaus);?
    
    hmass->Fit("fit2","W", "", 2.4, 4.);
    
    TF1 * expofit = new TF1("expofit","[0]*exp(-[1]*x)", 0., 1000.);//exponential function
    expofit->SetParameters(fit2->GetParameter(3), fit2->GetParameter(4));
    expofit->SetLineColor(kBlue);//Draw using different color
    
    cout << "*** fit to counts ***" << endl;
    cout << "expo bkg " << fit2->GetParameter(3) << "\t"
    // Draw backgrounds
    hmass->Draw();
    hmass->GetXaxis()->SetRangeUser(2.4,4.0);
    expofit->DrawCopy("same");
    //expo bkg + - 1 sigma error on bkg amplitude
    Double_t expobkg_sigmap = fit2->GetParameter(3) + fit2->GetParError(3);
    expofit->SetParameters(expobkg_sigmap, fit2->GetParameter(4));
    expofit->SetLineColor(kGreen);
    expofit->DrawCopy("same");
    
    Double_t expobkg_sigmam = fit2->GetParameter(3) - fit2->GetParError(3); //
    expofit->SetParameters(expobkg_sigmam, fit2->GetParameter(4));
    expofit->SetLineColor(kMagneta);
    expofit->Draw("same");
    
    return;
    
    cc->cd(2);
    int nbins = hmass->GetNbinsX();
    for (int ibin=1; ibin<=nbins; ibin++)
    {
        double x = hmass->GetXaxis()->GetBinCenter(ibin);
        double fval = expofit->Eval(x);
        hsub->SetBinContent(ibin,hmass->GetBinContent(ibin)-fval);
        
        cout << ibin << "\t hmass = " << hmass->GetBinContent(ibin) << "\tfval = " << fval  <<"\t hsub = "<< hsub->GetBinContent(ibin) << endl;
        
    }
    
    hsub->Draw("e");
    hsub->GetXaxis()->SetRangeUser(2.4,4.0);
    
    
    // +1 sigma
    cc->cd(3);
    int nbins = hmass->GetNbinsX();
    for (int ibin=1; ibin<=nbins; ibin++)
    {
        double x = hmass->GetXaxis()->GetBinCenter(ibin);
        double fval = expofit->Eval(x);
        hsub->SetBinContent(ibin,hmass->GetBinContent(ibin)-fval);
        cout << ibin << "\t hmass = " << hmass->GetBinContent(ibin) << "\tfval = " << fval  <<"\t hsub = "<< hsub->GetBinContent(ibin) << endl;
        
    }
    
    hsub->Draw("e");
    hsub->GetXaxis()->SetRangeUser(2.4,4.0);
    
    
    cc->cd(2);
    int nbins = hmass->GetNbinsX();
    for (int ibin=1; ibin<=nbins; ibin++)
    {
        double x = hmass->GetXaxis()->GetBinCenter(ibin);
        double fval = expofit->Eval(x);
        hsub->SetBinContent(ibin,hmass->GetBinContent(ibin)-fval);
    
        
        cout << ibin << "\t hmass = " << hmass->GetBinContent(ibin) << "\tfval = " << fval  <<"\t hsub = "<< hsub->GetBinContent(ibin) << endl;
        
    }
    
    hsub->Draw("e");
    hsub->GetXaxis()->SetRangeUser(2.4,4.0);
    
    
}