// #ifndef __CINT__
#include "RooGlobalFunc.h"
//#endif
#include "RooRealVar.h"
#include "RooArgList.h"
#include "RooFormulaVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooExponential.h"
#include "RooBifurGauss.h"
#include "RooAddModel.h"
#include "RooProdPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooHist.h"
#include "RooCBShape.h"
#include "RooPolynomial.h"
#include "RooBinning.h"
#include "TH1.h"
#include "TH2.h"
#include "RooAddPdf.h"
#include "RooProdPdf.h"
#include "RooFitResult.h"
#include "RooGenericPdf.h"
#include "RooLandau.h"
#include "TChain.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
#include "RooCategory.h"
#include "RooSuperCategory.h"
#include "RooSimultaneous.h"
#include "RooNLLVar.h"
#include "TFile.h"
#include "TPaveLabel.h"
#include "RooBinning.h"
#include "RooDataHist.h"
#include "RooJohnson.h"
#include "RooHypatia2.h"
#include "RooCrystalBall.h"

#include "/home/chanchal/Documents/chanchal phd/analysis_new/Dstatg_CFT_GenMC_feb_2023/fit/WS_final_fit/myRooJohnsonSU.cpp"
#include "/home/chanchal/Documents/chanchal phd/analysis_new/Dstatg_CFT_GenMC_feb_2023/fit/WS_final_fit/myRooPolBG.cpp"

using namespace RooFit ;
using namespace std;


#include "TLegend.h"
#include "TLatex.h"
#include "TROOT.h"
#include "RooWorkspace.h"
#include "TStyle.h"
#include "TGaxis.h"


void fit_d0pi_binsNew_twopdf_check()
{
 
    // // // // Open the ROOT file
    TFile* file = new TFile("RS_MC15ri_1ab_total_fit_BDT_binarypid_neweffcorr_pisoftc_400bind0pi.root");
    TH2 *h = (TH2*)file-> Get("2D"); 
   
    // // // // MD0 histogram
    TH1F* h_md0 = new TH1F("h_md0", "candidates that peak in Md0pi", 200, 1.8, 1.95);

    const int numBinsX = h->GetNbinsX();
    cout << " GetBinx = " << h->GetNbinsX() << " GetBiny = " << h->GetNbinsY() << endl;
       

    // // // Create an array of TH1D pointers
    TH1D* h_projections[numBinsX];


    RooRealVar d0pi("d0pi","", 2.0044, 2.022,"");


    // // // signal
    // /// // model 1
    RooRealVar f_G1("f_{G1}", "",0.68, 0, 1);  // 9, 6                                                                                             
    RooRealVar mu_G1("#mu_{G1}", "mean_gauss1", 2.01, 2.008, 2.012);                                                         
    RooRealVar sigma_G1("#sigma_{G1}", "sigma_gauss1", 0.05, 0, 0.1);                                                                       
    RooRealVar sigma_G2("#sigma_{G2}", "sigma_gauss2", 0.02, 0, 0.1);
    RooGaussian gauss1("gauss1", "1st gaussian PDF", d0pi, mu_G1, sigma_G1);                                                           
    RooGaussian gauss2("gauss2", "2nd gaussian PDF", d0pi, mu_G1, sigma_G2);                                                           
    RooAddPdf sig1("sig1", "double_gaussian", gauss1, gauss2, f_G1);


    
    // // model 2
    RooRealVar mu_J("#mu_{J}","mean_johnson", 2.01025, 2.0044, 2.022);
    RooRealVar sigma_J("#sigma_{J}", "sigma_johnson", 0.000273251, 0, 0.01);
    RooRealVar gamma_J("#gamma_{J}", "gamma", -0.0643552, -1, 1);
    RooRealVar delta_J("#delta_{J}", "delta", 1.06948, 0, 2);
    myRooJohnsonSU sig2("sig2", "johnson", d0pi, mu_J, sigma_J, delta_J, gamma_J);


  

    // // // Non peaking bkg
    RooRealVar threshold_comb("threshold_comb", "no data beyond this", 2.00441);
    RooRealVar alpha1_comb("#alpha1_comb", "alpha", -43.0458, -50, 50); 
    RooRealVar beta1_comb("#beta1_comb", "beta", 0);  
    myRooPolBG bkg1("bkg1","polynomial", d0pi, threshold_comb, alpha1_comb, beta1_comb);


    RooRealVar alpha2_comb("#alpha2_comb", "alpha", -49.0458, -80, 80); 
    RooRealVar beta2_comb("#beta2_comb", "beta", 0);  
    myRooPolBG bkg2("bkg2","polynomial", d0pi, threshold_comb, alpha2_comb, beta2_comb);


    // // // total
    RooRealVar n_sig1("N_{sig1}", "n_{s}", 11480, 0, 500000);
    RooRealVar n_bkg1("N_{bkg1}", "n_{b}", 9862, 0, 500000);

 
    RooRealVar n_sig2("N_{sig2}", "n_{s}", 175660, 0, 500000);
    RooRealVar n_bkg2("N_{bkg2}", "n_{b}", 17302, 0, 500000);

   
    // // // Define different models
    std::vector<RooAbsPdf*> models;
    models.push_back(new RooAddPdf("pdf1", "two component model 1", RooArgList(sig1, bkg1), RooArgList(n_sig1, n_bkg1)));
    models.push_back(new RooAddPdf("pdf2", "two component model 2", RooArgList(sig2, bkg2), RooArgList(n_sig2, n_bkg2)));
   

    // // // Create a TFile to save the histograms
    TFile* outputFile = new TFile("fit_results_RS.root", "RECREATE");

    // // Open a PDF file to store histograms
    TCanvas pdfCanvas("pdfCanvas", "Fitting Results", 800, 600);
    pdfCanvas.Print("fit_results_RS.pdf[");



    // // Loop over the bins along the x-axis
    for (int bin_x = 1; bin_x <= numBinsX; ++bin_x) 

    {
        
        cout << "///////////////////////////// bin no = " << bin_x << "///////////////////////////" << endl;
      

        // // // Project along the y-axis for the current bin along x-axis
        h_projections[bin_x - 1] = h->ProjectionY(Form("h1_projection_x%d", bin_x), bin_x, bin_x);
       
        RooDataHist data("data","data", RooArgList(d0pi), h_projections[bin_x - 1]);

      
        // // // Choose the model based on some criteria 
        int modelIndex = (bin_x > 80 && bin_x < 100) ? 1 : 0;  
        RooAbsPdf* pdf = models[modelIndex];

       // // Fit the model to the data
       RooFitResult *fitresult = pdf->fitTo(data, Save(true), Strategy(2), Extended(true), RooFit::SumW2Error(true));
       
       // Check the condition for setting n_sig_value
       double n_sig_value, n_sig_error, n_bkg_value, n_bkg_error;


        if (bin_x > 80 && bin_x < 100)
        {
          n_sig_value = n_sig2.getVal();
          n_sig_error = n_sig2.getError();
          
          n_bkg_value = n_bkg2.getVal();
          n_bkg_error = n_bkg2.getError();
        }
         else
        {
          n_sig_value = n_sig1.getVal();
          n_sig_error = n_sig1.getError();

          n_bkg_value = n_bkg1.getVal();
          n_bkg_error = n_bkg1.getError();
        }


       cout << "n_sig_value = " << n_sig_value << " +/- " << n_sig_error << endl;
       cout << "n_bkg_value = " << n_bkg_value << " +/- " << n_bkg_error << endl;

       
      // // // // // // Md0 that peak in Md0pi
      // h_md0 ->SetBinContent(bin_x, n_sig_value);
      // h_md0->SetBinError(bin_x, n_sig_error);


      h_md0 ->SetBinContent(bin_x, n_bkg_value);
      h_md0->SetBinError(bin_x, n_bkg_error);
     
     


      // // //  plot D0pi 
      RooPlot* d0piplot = d0pi.frame(Title(Form("D0pi distribution for Bin %d", bin_x)));
      data.plotOn(d0piplot);
      pdf->plotOn(d0piplot);
      pdf->paramOn(d0piplot);
      d0piplot->SetYTitle("Candidates per 0.09 MeV/#it{c}^{2}");
      d0piplot->GetYaxis()->CenterTitle();
      d0piplot->GetYaxis()->SetTitleSize(0.05);
      d0piplot->GetXaxis()->SetLabelOffset(999);
      // d0piplot->GetYaxis()->SetLogy(true);


      // // // chi prob
      cout << "chi^2 = " << d0piplot->chiSquare() << endl;
      double nbins = 11*d0piplot->chiSquare(11) / (d0piplot->chiSquare(11) - d0piplot->chiSquare());// 10 is the no. of fit parameters
      cout << "chi^2/ndf = " << d0piplot->chiSquare()*nbins << "/" << nbins-11 << endl;
      cout << "chi^2/ndf = " << d0piplot->chiSquare(11)<< endl;
      cout << "100bin_histo: chi^2/ndf = " << d0piplot->chiSquare()*100 << "/" << 100-11 << endl;
      
      double d0piplot_chi2 = d0piplot->chiSquare();
      TPaveLabel *t1 = new TPaveLabel(0.7,0.4,0.9,0.46, Form("#chi^{2}/ndf = %f", d0piplot_chi2),"brNDC");
      t1->SetBorderSize(1);
      t1->SetFillStyle(0);
      //t1->SetTextFont(40);
      t1->SetTextSize(0.70);
      t1->SetLineColor(0);


      // pull distribution
      RooHist* hpull = d0piplot->pullHist();
      hpull->SetFillStyle(1001);
      hpull->SetFillColor(1);
      for(int i=0;i<hpull->GetN();++i) hpull->SetPointError(i,0.0,0.0,0.0,0.0);
      RooPlot* pullplot = d0pi.frame(Title(" "));
      pullplot->addPlotable(hpull,"B");
      // pullplot->SetYTitle("Pull");
      pullplot->SetXTitle("m(D^{0}#pi^{+}) [GeV/#it{c}^{2}]");
      pullplot->SetMinimum(-4.);
      pullplot->SetMaximum(4.);
      pullplot->GetXaxis()->SetLabelSize(0.11);
      pullplot->GetXaxis()->SetTitleSize(0.13);
      pullplot->GetYaxis()->SetLabelSize(0.08);
      pullplot->GetXaxis()->CenterTitle();


      TCanvas *canvas = new TCanvas(Form("fitCanvas_x%d", bin_x), Form("Fit Result for Bin %d", bin_x), 800, 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);
      // canvas->cd(1)->SetLogy(true);
       
      pdf->plotOn(d0piplot,Components(RooArgSet(*pdf)),Name("pdf"));
      pdf->plotOn(d0piplot,Components(RooArgSet(sig1, sig2)),LineColor(kRed),LineStyle(kDashed), Name("sig"));
      pdf->plotOn(d0piplot,Components(RooArgSet(bkg1, bkg2)),LineColor(kBlack),LineStyle(kDashed), Name("bkg"));
   
   
      d0piplot->Draw();
      t1->Draw();
      // upPad->SetLogy(1);

      canvas->cd(2);
      pullplot->Draw();
      
      
      // canvas->Update();
   
     // // Save the canvas as a PDF
    canvas->Print("fit_results_RS.pdf");

    // Save the canvas to the output file
    outputFile->cd();
    canvas->Write();

    delete canvas;

     }
   

    // // // Close the PDF file and TFile
    pdfCanvas.Print("fit_results_RS.pdf]");
    outputFile->Close();
       
 
     // h_md0 ->Draw();

     

  TFile * f_output = new TFile("Md0_genericMC_peakinmd0pi_twomodels_400bins_binerror.root","recreate");
  h_md0->Write(); 
  f_output->Close();  
      
   
}