#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "TFile.h"
#include <cmath>
 
double hA1( double w , double rho2, double r11, double r21);
double R1w( double w , double rho2, double r11, double r21);
double R2w( double w );

double h_p(double w, double rho2, double r11, double r21);
double h_m(double w, double rho2, double r11, double r21);
double h_0(double w, double rho2, double r11, double r21);

double w_diff1(double w, double rho2, double r11, double r21);
double w_integration1(double min, double max, double rho2, double r11, double r21);

double mB  = 5.27958;
double mD  = 2.01022;
double r = mD/mB;
double rr = (2*sqrt(mB*mD)) /(mB+mD);
double rho2 = 1.3;
double r11 = 1.18;
double r21 = 0.71;
double hA11 =0.921 ;
double Vcb = 0.0414; 
double gF = 1.16637E-5;
double factor = (( 2 * TMath::Pi()) * (((gF*gF)*(Vcb*Vcb) * (6 * mB * mD *mD))  / (2048 * (TMath::Pi()*TMath::Pi()*TMath::Pi()* TMath::Pi()) )));
double constant =0;
int nbins = 10;
int N = 30;

void ratiotest()
{

	TFile* f3   = new TFile("/home/waheede/phd_analysis/test_mc/mc_test.root");//signal MC

	double constant = ( (152e+6 ) * 0.677 * 0.0382 * 1.519e-12 ) / (6.58212E-25); //SVD1

	double minw, maxw;
	minw = 1.0;  maxw = 1.0;

	//definition of input histograms
	TH1F* h_w1_up = new TH1F("w1_up", "w; w; Events", 10,1.0,1.504);
	TH1F *hw;
        f3->GetObject("w" , hw);
	
	//ratio histograms for each method i tried
	TH1F * h_w_r1;
	TH1F * h_w_r12;
	TH1F * h_w_r13;
        h_w_r1 = new TH1F("w_r1", "w; w ; Ratio", 10,1,1.504);
	h_w_r12 = new TH1F("w_r12", "w; w ; Ratio", 10,1,1.504);
	h_w_r13 = new TH1F("w_r13", "w; w ; Ratio", 10,1,1.504);

        for(int j = 0; j<10; j++){
                maxw      = maxw + 0.054;
		h_w1_up->SetBinContent(j+1,  constant * abs(w_integration1(minw, maxw, rho2, r11, r21)));
                minw      = maxw;
             }


	TCanvas *c = new TCanvas("c","Efficiency",900,400);
        c->Print("ratiotest.pdf[");
	
        TLine *line2 =new TLine(1,1,1.504,1);
        line2->SetLineColor(kRed);

	//Hisotgrams with Normalization
	h_w1_up->SetLineColor(kRed);
	hw->SetLineColor(kBlue);
   	h_w1_up->Scale(1./h_w1_up->Integral()); 
   	hw->Scale(1./hw->Integral()); 
 	hw->SetMinimum(0);
	h_w1_up->SetMinimum(0);
   	hw->SetStats(0);
   	h_w1_up->SetStats(0);
	h_w1_up->Draw("hist");
        hw->Draw("same");	
	c->Print("ratiotest.pdf",  "Title: "  "w Theory ");

	//Ratio with Divide(drawing with hist)
	 h_w1_up->Scale(1./h_w1_up->Integral());
        hw->Scale(1./hw->Integral());
	h_w_r1->Divide(hw,h_w1_up, 1,1);
	h_w1_up->SetMinimum(0.9);
	h_w1_up->SetMaximum(1.1);	
	h_w_r1->Draw("hist");
	line2->Draw();
	c->Print("ratiotest.pdf",  "Title: "  "w Theory ");

	//Ratio with Divide(drawing with error option)
	h_w1_up->Scale(1./h_w1_up->Integral());
        hw->Scale(1./hw->Integral());
	h_w_r1->Divide(hw,h_w1_up, 1,1);
 	h_w1_up->SetMinimum(0.9);
        h_w1_up->SetMaximum(1.1);
	h_w_r1->SetMarkerStyle(21);
        h_w_r1->Draw("e");
	line2->Draw();
        c->Print("ratiotest.pdf",  "Title: "  "w Theory ");

	//Using TGraph
	hw->Sumw2();
   	h_w1_up->Sumw2();
	h_w1_up->Scale(1./h_w1_up->Integral());
        hw->Scale(1./hw->Integral());

	auto ratio = new TGraphAsymmErrors(50); 
   	ratio->Divide(hw,h_w1_up,"pois"); 
   	ratio->SetMarkerStyle(20); 
   	ratio->GetXaxis()->SetRangeUser(1,1.504); 
   	ratio->SetMinimum(0.9);
   	ratio->SetMaximum(1.1); 
   	ratio->Draw("AP");
   	line2->Draw();
	c->Print("ratiotest.pdf",  "Title: "  "w Theory ");

   	ratio->SetTitle("Ratio using Poisson statistics");
   	ratio->Draw("P");
	c->Print("ratiotest.pdf",  "Title: "  "w Theory ");

	//trying to add the error
//	for (int i=1; i<=hw->GetNbinsX(); i++)
  //      {
//	ratio->SetPointError(i, 0,0,(h_w1_up->GetBinContent(i)/hw->GetBinContent(i)) );
//	cout<"error:"<<();
  //      }

	float red_histo_integral = hw->Integral();
	float blue_histo_integral = h_w1_up->Integral();
	float rel_uncert_overall = 0;
	for (int i=1; i<=(h_w_r12->GetNbinsX()); i++)
	{
    	float rel_uncert_from_blue_histo = (hw->GetBinContent(i) == 0)?0.:hw->GetBinError(i)/blue_histo_integral/hw->GetBinContent(i);
    	float rel_uncert_from_red_histo = (h_w1_up->GetBinContent(i) == 0)?0.:h_w1_up->GetBinError(i)/red_histo_integral/h_w1_up->GetBinContent(i);
    	rel_uncert_overall = sqrt(pow(rel_uncert_from_blue_histo,2)+pow(rel_uncert_from_red_histo,2));
    	h_w_r12->SetBinError(i, rel_uncert_overall*h_w_r12->GetBinContent(i));
	}
	h_w_r12->Draw("e");
	c->Print("ratiotest.pdf",  "Title: "  "w Theory ");


	//looking at stat. errors of both histograms (again, bin-by-bin)
   	for (i=1; i<=hw->GetNbinsX(); i++)
      	printf("Stat. errors for non-scaled histograms: bin #%02i: red = %07.3f%%, blue = %07.3f%%\n", i, h_w1_up->GetBinError(i)*100/h_w1_up->GetBinContent(i), hw->GetBinError(i)*100/hw->GetBinContent(i));

   	h_w1_up->SetMaximum(h_w1_up->GetMaximum()*1.5);

	//getting integrals before scaling histograms:
   	hw->Sumw2();
   	h_w1_up->Sumw2();

   	printf("Now scaling histograms...\n");

   	hw->Scale(1/hw->Integral());
   	h_w1_up->Scale(1/h_w1_up->Integral());

   	for (int i=1; i<=hw->GetNbinsX(); i++)
   	{
      	printf("Stat. errors for scaled histograms: bin #%02i: red = %07.3f%%, blue = %07.3f%%\n", i, h_w1_up->GetBinError(i)*100/h_w1_up->GetBinContent(i), hw->GetBinError(i)*100/hw->GetBinContent(i));
   	}

   	h_w1_up->Draw("hist");
   	hw->Draw("same");
	c->Print("ratiotest.pdf",  "Title: "  "w Theory ");

	//now dividing histograms:
	h_w_r13->Divide(hw,h_w_r1);
   	h_w_r13->GetYaxis()->SetTitleOffset(1.2);
   	h_w_r13->SetStats(0);

   	h_w_r13->SetMarkerStyle(21);
   	h_w_r13->SetMaximum(1.1);
   	h_w_r13->SetMinimum(0.9);
   	h_w_r13->Draw("ep");
	line2->Draw();
 	c->Print("ratiotest.pdf",  "Title: "  "w Theory ");

//now looking at the stat. errors of this ratio histogram. They equal sqrt(rel_unc_of_histo_in_numerator^2 + rel_unc_of_histo_in_denominator^2)

   	for (int i=1; i<=h_w_r1->GetNbinsX(); i++)
      	printf("Stat. errors for a ratio histograms: bin #%02i: %07.3f%%\n", i, h_w_r1->GetBinError(i)*100/h_w_r1->GetBinContent(i));

	c->Print("ratiotest.pdf",  "Title: "  "w Theory ");	
        c->Print("ratiotest.pdf]");

}

double w_integration1(double min, double max, double rho2, double r11,double r21)
{
double ans = 0;
int n = 21;
double h = (max - min)/n;
double sum1 = w_diff1( min + (h/2), rho2,r11, r21);
double sum2 = 0;
for (int i = 1; i<(n-1); i++ )
        {
                        sum1 = sum1 + w_diff1((min + (h * i) + h/2), rho2,r11, r21);
                        sum2 = sum2 + w_diff1((min + (h * i)), rho2,r11, r21);
        }
ans = (h / 6) * ( w_diff1(min, rho2,r11, r21)+ w_diff1(max, rho2,r11, r21)+ (4 *sum1) +( 2 * sum2));
return ans;
}

double w_diff1(double w,  double rho2, double r11,double r21)
{

double result = 0 ;
double f1 = (((gF*gF)*(Vcb*Vcb) * (mD * mD *mD)* ((mB-mD)*(mB-mD)))/(48 * TMath::Pi()*TMath::Pi()*TMath::Pi()));
double f2 = (sqrt((w*w)-1)) * ((w +1)*(w+1)) * (hA1(w,rho2, r11, r21) * hA1(w,rho2, r11, r21));
double f3 = h_p(w,rho2,r11,r21) + h_m(w,rho2,r11,r21) + h_0(w,rho2,r11,r21);
result = f1 *  f2 * f3 ;
return abs(result);
}

double z( double w )
{
double ans=0;
ans = (sqrt(w+1)-sqrt(2)) / ( sqrt(w+1)+sqrt(2));
return ans;
}
double R1w( double w , double rho2, double r11,double r21)
{
return ( r11 - ((0.12)* (w - 1)) + (0.05 * (w - 1) * ( w - 1)) );
}
double R2w( double w , double rho2, double r11,double r21)
{
return ( r21 + ((0.11) * (w - 1)) - (0.06 * (w - 1) * ( w - 1) ));
}
double hA1( double w , double rho2, double r11,double r21)
{
double ans = 0;
ans = hA11 * ( 1 -  ( 8 * rho2 * z(w)) + (( 53 * rho2 - 15 )*(z(w) * z(w))) -  (((231 * rho2) - 91) * (z(w) * z(w) * z(w))));
return ans;
}
double h_p(double w, double rho2, double r11,double r21)
{
double ans = 0.;
ans =   ((1 -  (2 * w * r) + (r * r))  / (( 1 - r) * (1-r)))  * (( 1 -  (sqrt ( (w - 1)/(w + 1) ) * R1w (w,rho2,r11,r21)) ) * ( 1 -  (sqrt ( (w - 1)/(w + 1) ) * R1w (w,rho2,r11,r21)) ) );
return ans;
}
double h_m(double w, double rho2, double r11,double r21)
{
double ans = 0.;
ans =  ( (1 -  (2 * w * r) + (r * r))  / (( 1 - r) * (1-r)))  * (( 1 +  (sqrt ( (w - 1)/(w + 1) ) * R1w (w,rho2,r11,r21)) ) * ( 1 +  (sqrt ( (w - 1)/(w + 1) ) * R1w (w,rho2,r11,r21)) ) );
return ans;
}
double h_0(double w,double rho2, double r11, double r21)
{
double ans = 0.;
ans = (1 + ( (w - 1)/(1-r) ) * ( 1 - R2w ( w ,rho2,r11,r21)) ) * ( 1 + (( (w - 1)/(1-r) ) * ( 1 - R2w ( w ,rho2,r11,r21)) ));
return ans;
}