#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TPad.h>
#include "TLorentzVector.h"
#include <THistPainter.h>
#include "THStack.h"
#include"TFile.h"
#include "TH1.h"
// put here the include files you need                                                                                                                                             
#include "TStopwatch.h"
#include <iostream>
using namespace std;
#include <iomanip>
#include <fstream>
#include <string>
#include <sstream>
using namespace std;
void Stack_Hist_Signals()

{
    double lumi = 300;
    gStyle->SetOptStat(0); //to switch for stat. box
  TString dirname = "/home/ahmed/bb-quark_Signals";
  TString rootH0BB = "h90BP1_bb.root";
  TString rootH1BB= "h90BP2_bb.root";
  TString rootH2BB  = "h95BP1_bb.root";
  TString rootH3BB  = "h95BP3_bb.root";
  TString rootH4BB  = "h96BP1_bb.root";
  TString rootH5BB  = "h96BP2_bb.root";
  TString rootH6BB  = "h96BP3_bb.root";
  TString rootH7BB = "h96BP6_bb.root";
  TString rootH8BB = "h_bb.root";
  TString rootH9BB = "z_bb.root";
  TString rootH10BB = "bb.root";            
  //TString ps_name = "bb_plots_log.ps";    
 
  TString ps_name = "bb_plots_log.ps";   

   rootH0BB=dirname + "/" + rootH0BB;
  
   rootH1BB=dirname + "/" + rootH1BB;
   rootH2BB=dirname + "/" + rootH2BB;
   rootH3BB =dirname + "/" + rootH3BB ;
   rootH4BB =dirname + "/" + rootH4BB ;
   rootH5BB =dirname + "/" + rootH5BB ;
   rootH6BB =dirname + "/" + rootH6BB ;
   rootH7BB =dirname + "/" + rootH7BB ;
   rootH8BB =dirname + "/" + rootH8BB ;
   rootH9BB =dirname + "/" + rootH9BB ;   
   rootH10BB =dirname + "/" + rootH10BB ;     
   

  delete gROOT->GetListOfFiles()->FindObject(rootH0BB); 
 delete gROOT->GetListOfFiles()->FindObject( rootH1BB); 
  delete gROOT->GetListOfFiles()->FindObject(rootH2BB); 
  delete gROOT->GetListOfFiles()->FindObject(rootH3BB);
  delete gROOT->GetListOfFiles()->FindObject(rootH4BB);
  delete gROOT->GetListOfFiles()->FindObject(rootH5BB);
  delete gROOT->GetListOfFiles()->FindObject(rootH6BB);
  delete gROOT->GetListOfFiles()->FindObject(rootH7BB);
  delete gROOT->GetListOfFiles()->FindObject(rootH8BB);
  delete gROOT->GetListOfFiles()->FindObject(rootH9BB);  
 
 
  TFile *rootH0BBbar= new TFile(rootH0BB);
 TFile *rootH1BBbar = new TFile(rootH1BB);
  TFile *rootH2BBbar= new TFile(rootH2BB);
  TFile *rootH3BBbar = new TFile(rootH3BB);
  TFile *rootH4BBbar = new TFile(rootH4BB);
  TFile *rootH5BBbar = new TFile(rootH5BB);
  TFile *rootH6BBbar = new TFile(rootH6BB);
  TFile *rootH7BBbar = new TFile(rootH7BB);
  TFile *rootH8BBbar = new TFile(rootH8BB);
  TFile *rootH9BBbar = new TFile(rootH9BB);    
  TFile *rootH10BBbar = new TFile(rootH10BB); 
                                              
                                               
  
  TH1F *hmc[11][9];TString hname_mc[9];
  
 hname_mc[0]="b1_pt";
 hname_mc[1]="b1_eta";
  hname_mc[2]="b2_pt";
  hname_mc[3]="b2_eta";
  hname_mc[4]="h_mass";
  hname_mc[5]="h_pt";
  hname_mc[6]="h_eta";
  hname_mc[7]="DR_b";  
  hname_mc[8]="scalarHt_Ht";
 
     
     
           


   
 
  
   
  
 
  
   
 //double Y_min[8]; double Y_max[8];

       for(int i=0; i<9;i++ ){

    //cout <<"... coloring mc hmc["<<i<<"] " << endl;
    
     hmc[0][i] = (TH1F*)rootH0BBbar->Get(hname_mc[i]);
    hmc[1][i] = (TH1F*)rootH1BBbar->Get(hname_mc[i]);
     hmc[2][i] = (TH1F*)rootH2BBbar->Get(hname_mc[i]);
     hmc[3][i] = (TH1F*)rootH3BBbar->Get(hname_mc[i]);
     hmc[4][i] = (TH1F*)rootH4BBbar->Get(hname_mc[i]);
     hmc[5][i] = (TH1F*)rootH5BBbar->Get(hname_mc[i]);
     hmc[6][i] = (TH1F*)rootH6BBbar->Get(hname_mc[i]);
     hmc[7][i] = (TH1F*)rootH7BBbar->Get(hname_mc[i]);
     hmc[8][i] = (TH1F*)rootH8BBbar->Get(hname_mc[i]);
     hmc[9][i] = (TH1F*)rootH9BBbar->Get(hname_mc[i]); 
      hmc[10][i] = (TH1F*)rootH10BBbar->Get(hname_mc[i]);        
  
/*
 Y_min[i] =(hmc[0][i]->GetMaximum()/1000);
    if(i==1||i==3||i==6)Y_max[i] = 100000*(hmc[2][i]->GetMaximum());
   else if(i==2||i==4||i==7)Y_max[i] = 100000*(hmc[2][i]->GetMaximum());
  else if(i==6)Y_max[i] = 100000*(hmc[2][i]->GetMaximum());
    else Y_max[i] = 100000*(hmc[2][i]->GetMaximum());

    hmc[0][i]->GetYaxis()->SetRangeUser(Y_min[i],Y_max[i]);
     if(i==0||i==2||i==4|| i==5)hmc[0][i]->GetXaxis()->SetRangeUser(0.,400.);
     if(i==1||i==3|| i==6||i==7)hmc[0][i]->GetXaxis()->SetRangeUser(-9.,9.);
     if(i==0||i==2||i==4|| i==5)hmc[1][i]->GetXaxis()->SetRangeUser(0.,400.);
     if(i==1||i==3|| i==6||i==7)hmc[1][i]->GetXaxis()->SetRangeUser(-9.,9.);
     if(i==0||i==2||i==4|| i==5)hmc[2][i]->GetXaxis()->SetRangeUser(0.,400.);
     if(i==1||i==3|| i==6||i==7)hmc[2][i]->GetXaxis()->SetRangeUser(-9.,9.);
     if(i==0||i==2||i==4|| i==5)hmc[3][i]->GetXaxis()->SetRangeUser(0.,400.);
     if(i==1||i==3||i==7)hmc[3][i]->GetXaxis()->SetRangeUser(-9.,9.);
  
  */

     hmc[0][i]->SetLineColor(kYellow);
   //  hmc[0][i]->SetFillColor(kYellow );

    
      hmc[1][i]->SetLineColor(kViolet );
    //  hmc[1][i]->SetFillColor(kViolet);
      
      hmc[2][i]->SetLineColor(kBlue);
     // hmc[2][i]->SetFillColor(kBlue);
      
     hmc[3][i]->SetLineColor(kGray);
    // hmc[3][i]->SetFillColor(kGray);
  
  
     hmc[4][i]->SetLineColor(kRed -9 );
   //  hmc[4][i]->SetFillColor(kRed -9);
    
     hmc[5][i]->SetLineColor(kRed -6);
   //  hmc[5][i]->SetFillColor(kRed -6);
     
     hmc[6][i]->SetLineColor(kRed -10);
   //  hmc[6][i]->SetFillColor(kRed-10);
     
     hmc[7][i]->SetLineColor(kRed );
   //  hmc[7][i]->SetFillColor(kRed );
     
      hmc[8][i]->SetLineColor(kOrange );
   //  hmc[8][i]->SetFillColor(kOrange );
   
    hmc[9][i]->SetLineColor(kCyan );
   //  hmc[9][i]->SetFillColor(kCyan );   
   
    hmc[10][i]->SetLineColor(kGreen );
   //  hmc[10][i]->SetFillColor(kGreen );
    
     
 
  
     hmc[0][i]->SetMarkerColor(kYellow);
    
    hmc[1][i]->SetMarkerColor( kViolet);  
     hmc[2][i]->SetMarkerColor(kBlue);  
     hmc[3][i]->SetMarkerColor(kGray);
     hmc[4][i]->SetMarkerColor(kRed -9);
     hmc[5][i]->SetMarkerColor(kRed -6);
     hmc[6][i]->SetMarkerColor(kRed -10);
     hmc[7][i]->SetMarkerColor(kRed );
     hmc[8][i]->SetMarkerColor(kOrange );
     hmc[9][i]->SetMarkerColor(kCyan);    
     hmc[10][i]->SetMarkerColor(kGreen);      
    
   
     hmc[0][i]->Scale(1/hmc[0][i]->Integral());
     hmc[1][i]->Scale(1/hmc[1][i]->Integral());
     hmc[2][i]->Scale(1/hmc[2][i]->Integral()); 
     hmc[3][i]->Scale(1/hmc[3][i]->Integral());
     hmc[4][i]->Scale(1/hmc[4][i]->Integral());
     hmc[5][i]->Scale(1/hmc[5][i]->Integral());
     hmc[6][i]->Scale(1/hmc[6][i]->Integral());
     hmc[7][i]->Scale(1/hmc[7][i]->Integral());
     hmc[8][i]->Scale(1/hmc[8][i]->Integral());
     hmc[9][i]->Scale(1/hmc[9][i]->Integral());  
     hmc[10][i]->Scale(1/hmc[10][i]->Integral());   
  
  
    }//mc for
      
    
   // cout<<" ...1 "<<endl;
    
  
   
    // THStack* t[8];
  
  //for(int j=0;j<8;j++)
    { 
    // t[j] = new THStack( "Stack_Hist",""); 
                                             
    
   //  printf("0 %d %g\n",j,hmc[0][j]->GetMaximum());
     //   t[j]->Add(hmc[0][j]);
    
    
   //   printf("1 %d %g\n",j,hmc[1][j]->GetMaximum());
     // t[j]->Add(hmc[1][j]);
         
     // printf("2 %d %g\n",j,hmc[2][j]->GetMaximum());
     // t[j]->Add(hmc[2][j]);
     // printf("\n");
}
     
    
     
  
  
 // cout<<" ...2 "<<endl;
  //TLegend* legend = new TLegend(0.78,0.835,0.98,0.995);
  
    TLegend* legend = new TLegend(0.7,0.55,0.88,.88,Form("L = %1.1lf fb^{-1}",lumi));

    legend->AddEntry( hmc[0][0],"h90BP1_bb ","signal0");//H1BBbar sample
    legend->AddEntry( hmc[1][0],"h90BP2_bb","signal1");//H2BBbar sample
    legend->AddEntry( hmc[2][0],"h95BP1_bb","signal2");//H3BBbar sample
    legend->AddEntry( hmc[3][0],"h95BP3_bb","signal3");//H4BBbar sample
    legend->AddEntry( hmc[4][0],"h96BP1-bb","signal4");//H5BBbar sample
    legend->AddEntry( hmc[5][0],"h96BP2_bb","signal5");//H6BBbar sample
    legend->AddEntry( hmc[6][0],"h96BP3_bb","signal6");//H7BBbar sample
    legend->AddEntry( hmc[7][0],"h96BP6_bb","signal7");//H8BBbar sample
    legend->AddEntry( hmc[8][0],"h_bb","f0");          //H9BBbar sample
    legend->AddEntry( hmc[9][0],"z_bb","f1");         //H10BBbar sample  
    legend->AddEntry( hmc[10][0],"pp_bb","f2");       //H11BBbar sample         
 
 
 
  legend->SetBorderSize(1.); //no border for legend 
  legend->SetFillColor(0);
  legend->SetHeader("L_{data} = 300 fb^{-1}");

 TPostScript* ps = new TPostScript(ps_name,112);//111 portrait 112 landscape 113 eps
  ps->Range(20.0,20.0);

  TCanvas* C = new TCanvas("Fig","Fig");
  gPad->SetLogy();

 for(int i=0; i<9; i++ ){
   //  if(i!=4)continue;

    //else gPad->SetLogy();

   
    //gPad->SetLogy(1);
    // t[i]->Draw("Hist,same");
      
  
            hmc[0][i]->Draw("Hist");
            hmc[1][i]->Draw("Hist,same") ;
            hmc[2][i]->Draw("Hist,same") ;
            hmc[3][i]->Draw("Hist,same") ;
            hmc[4][i]->Draw("Hist,same") ;
            hmc[5][i]->Draw("Hist,same") ;
            hmc[6][i]->Draw("Hist,same") ;
            hmc[7][i]->Draw("Hist,same") ;
            hmc[8][i]->Draw("Hist,same") ;
            hmc[9][i]->Draw("Hist,same") ;
            hmc[10][i]->Draw("Hist,same") ;  
     
     gPad->RedrawAxis();
     /* if(i==174)leg->Draw();
	else*/ legend->Draw();
     C->Update(); 
 }

 ps->Print();
 ps->Close();
 system("ps2pdf bb_plots_log.ps");cout<<" ps converted..."<<endl;
 /////////////////////////////////////////////////////////////////
 
 
double Ns1=hmc[0][4]->Integral();
double Ns2=hmc[1][4]->Integral();
double Ns3=hmc[2][4]->Integral();
double Ns4=hmc[3][4]->Integral();
double Ns5=hmc[4][4]->Integral();
double Ns6=hmc[5][4]->Integral();
double Ns7=hmc[6][4]->Integral();
double Ns8=hmc[7][4]->Integral();

double Nb1=hmc[8][4]->Integral();
double Nb2=hmc[9][4]->Integral();
double Nb3=hmc[10][4]->Integral();

cout<<" number of signal events h1mass 90BP1 = "<< Ns1<<endl;
cout<<" number of signal events h1mass 90BP2 = "<< Ns2<<endl;
cout<<" number of signal events h1mass 95BP1 = "<< Ns3<<endl;
cout<<" number of signal events h1mass 95BP3 = "<< Ns4<<endl;
cout<<" number of signal events h1mass 96BP1 = "<< Ns5<<endl;
cout<<" number of signal events h1mass 96BP2 = "<< Ns6<<endl;
cout<<" number of signal events h1mass 96BP3 = "<< Ns7<<endl;
cout<<" number of signal events h1mass 96BP6 = "<< Ns8<<endl;
cout<<"+++++++++++++++++++++++++++++++++++++++++"<<endl;
cout<<"+++++++++++++++++++++++++++++++++++++++++"<<endl;
cout<<" number of BG events H2 = " << Nb1<<endl;
cout<<" number of BG events Z  = " << Nb2<<endl;
cout<<" number of BG events BB = " << Nb3<<endl;
cout<<"+++++++++++++++++++++++++++++++++++++++++"<<endl;



cout<<" Sum of all BG = "<<Nb1+Nb2+Nb3<<endl;
cout<<"+++++++++++++++++++++++++++++++++++++++++"<<endl;

cout<<"+++++++++++++++++++++++++++++++++++++++++"<<endl;
double S1= Ns1/sqrt(Nb1+Nb2+Nb3);
cout<<" Signal 90BP1 signficance = " << S1<<endl;

double S2= Ns2/sqrt(Nb1+Nb2+Nb3);
cout<<" Signal 90BP2 signficance = " << S2<<endl;

double S3= Ns3/sqrt(Nb1+Nb2+Nb3);
cout<<" Signal 95BP1 signficance = " << S3<<endl;

double S4= Ns4/sqrt(Nb1+Nb2+Nb3);
cout<<" Signal 95BP3 signficance = " << S4<<endl;

double S5= Ns5/sqrt(Nb1+Nb2+Nb3);
cout<<" Signal 96BP1 signficance = " << S5<<endl;

double S6= Ns6/sqrt(Nb1+Nb2+Nb3);
cout<<" Signal 96BP2 signficance = " << S6<<endl;

double S7= Ns7/sqrt(Nb1+Nb2+Nb3);
cout<<" Signal 96BP3 signficance = " << S7<<endl;

double S8= Ns8/sqrt(Nb1+Nb2+Nb3);
cout<<" Signal 96BP6 signficance = " << S8<<endl;

}