#if !(defined(__CINT__) || defined(__CLING__)) || defined(__ACLIC__)
#include <iostream>
#include <fstream>
#include <sstream>
#include <TFile.h>
#include <TDirectory.h>
#include <TH2.h>
#include <TCanvas.h>
#include <TLorentzVector.h>
#include <TMath.h>
#include <TLegend.h>
#include <TLatex.h>
#include <TROOT.h>
#include "TROOT.h"
#include "TStyle.h"
#include "TDirectory.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TAxis.h"
#include "TPaveLabel.h"
#include "TPaveText.h"
#include "TText.h"
#include "TLorentzVector.h"
#include "TCanvas.h"
#include "TROOT.h"
#include "TAxis.h"
#include "TH1.h"
#include "TTree.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH1I.h"
#include "TCanvas.h"
#include "TLine.h"
#include "TMath.h"
#include "TVector3.h"
#include "TString.h"
#include "TLegend.h"
#include "TStyle.h"
#include "TSystem.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TH1.h>
#include <TF1.h>
#include <TLegend.h>
#include <TMath.h>
#include <TPad.h>
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include "TStyle.h"
#include "TPaveStats.h"
#include "TText.h"
#include "THStack.h"
#include "TCut.h"
#include "TAxis.h"
#include "TChain.h"
#include "tnp_muon_UPC_PbPb2018.h"
#include "TUnfoldDensity.h"
#endif
#include "ReadTauTauTree.C"
using namespace std;
const int nSample = 1;
//const char *sample[nSample]={"Data","QED_SC","QED_SL"};
const char *Sample[nSample]={"TauTau"};

void make_canvas(TCanvas *&);
void make_canvas_ratio(TCanvas *&);
void make_hist(TH1D *&, Color_t , int );
void make_hist_ratio(TH1D *&,Color_t, int) ;

const double luminosity       = 1642.9;

const double nEventsGeneratedGUPC = 5270000; // older number with less files 67262800 superchic 
const double xsecGeneratedGUPC    = 1060.8; // ub
double scaleFactors = 0.85 *  // NEE    21.12.2021
                        0.952; // electron reco+ID 21.12.2021


const double nEventsGeneratedSC = 1000000;
const double xsecGeneratedSC = 570;

double lumiNormGUPC = xsecGeneratedGUPC*luminosity*scaleFactors/nEventsGeneratedGUPC;
double lumiNormSC = xsecGeneratedSC*luminosity*scaleFactors/nEventsGeneratedSC;
//double norm_exgg = 1.09*xsecGeneratedLbL*luminosity/nEventsGeneratedLbL;  //31.12.21 (SF = 1.048, took the square here). 

void printOutput(TH1D *hist);

const double wt[nSample] = {lumiNormGUPC};
double W = 1.0;
double W2 = 1.0;
float getAcoBinSF(float acop);


void drawText(const char *text, float xp, float yp, int size){
  TLatex *tex = new TLatex(xp,yp,text);
  tex->SetTextFont(43);
  tex->SetTextSize(size);
  tex->SetTextColor(kMagenta);
  tex->SetLineWidth(1);
  tex->SetNDC();
  tex->Draw();
}
const int nMassbins=4;
double Massbin[nMassbins]={5.0,8.0,11.0,16.0};
const int nMassbin= sizeof(Massbin)/sizeof(double) - 1;


TH1D* subtractBkg(TH1D *hdata, TH1D *hmumufsr, const char* name, int nbins, float xmin, float xmax);

TH1D* getXSecHist(TH1D *hist, const char* name, int nbins, float xmin, float xmax, float lumi);

TCanvas* PlotHistsAndRatio(TCanvas* , TH1D* , TH1D*, TH1D*, double , double , double , double , double , double , const char *, bool); 
void makeUnfold_data(bool GenScale, bool RecoScale){
 
  //gROOT->LoadMacro("CMS_lumi.C");
  bool outOfFrame    = false;
  TH1::SetDefaultSumw2();
  gStyle->SetOptStat(0);
  


  TFile *outf= new TFile(Form("unfolding_histograms_Data.root"),"recreate");  
  
  TChain *TauTau[nSample];
  
  TH1D* hMu_Pt[nSample], *hEleMu_Mass[nSample], *hEleMu_vSumPt[nSample], *hAcoplanarity[nSample]; TH1F *hSF[nSample];

  TH1D* hTau_GenPt[nSample], *hMu_Gen_Pt[nSample], *hEleMu_Gen_Mass[nSample], *hEleMu_Gen_vSumPt[nSample];
  TH2D* hEleMu_RecoGenPt[nSample], *hdipho_RecoGen[nSample], *hEleMu_RecoGenInvmass[nSample], *hEleMu_RecoGenvSumPt[nSample],  *hEleMu_RecoGenTauPt[nSample];

  
  cout << "define hists"<< endl;
  
  int i =0;
  //for (int i = 0; i < nSample; i++){
    TauTau[i] = new TChain("output_tree");
//     TauTau[i]->Add("RootFile/mc_GammaUPC_EleMu_ArashNew_Test2_merged.root");
//    TauTau[i]->Add("RootFile/mc_GammaUPC_EleMu_ArashNew_20thSep_merged.root");
//    TauTau[i]->Add("RootFile/GenMatch/gUPC_Tree.root");
//    TauTau[i]->Add("RootFile/GenMatch/SC_Arash_Tree.root");
    TauTau[i]->Add("RootFile/GenMatch/mc_GammaUPC_TUnfoldTree_7thOctober_hasonlyOneMuonRecoAndGen_test_merged.root");
    hMu_Pt[i]    = new TH1D(Form("hMu_Pt%s", Sample[i]),"",12,2.5,14.5);
    //hMu_Pt[i]    = new TH1D(Form("hMu_Pt%s", Sample[i]),"",nPtbin,Ptbin);
    hEleMu_Mass[i]   = new TH1D(Form("hEleMu_Mass%s",Sample[i]),"",5,4,24);
   // hEleMu_Mass[i]   = new TH1D(Form("hEleMu_Mass%s",Sample[i]),"",nMassbin,Massbin);
    hEleMu_vSumPt[i]   = new TH1D(Form("hEleMu_vSumPt%s",Sample[i]),"",5,0,10);
    hAcoplanarity[i] = new TH1D(Form("hAcoplanarity%s",Sample[i]), "",40,0, 0.4);

    hTau_GenPt[i]         = new TH1D(Form("hTau_GenPt%s", Sample[i]),"",6,2.5,14.5);
  //hTau_GenPt[i]         = new TH1D(Form("hTau_GenPt%s", Sample[i]),"",nPtbin,Ptbin);
    hMu_Gen_Pt[i]   = new TH1D(Form("hMu_Gen_Pt%s",Sample[i]),"",6,2.5,14.5);
//  hMu_Gen_Pt[i]   = new TH1D(Form("hMu_Gen_Pt%s",Sample[i]),"",nPtbin,Ptbin);
    hEleMu_Gen_Mass[i]          = new TH1D(Form("hEleMu_Gen_Mass%s",Sample[i]),"",5,4,24);
    //hEleMu_Gen_Mass[i]          = new TH1D(Form("hEleMu_Gen_Mass%s",Sample[i]),"",nMassbin,Massbin);
    hEleMu_Gen_vSumPt[i]          = new TH1D(Form("hEleMu_Gen_vSumPt%s",Sample[i]),"",5,0,10);
    hEleMu_RecoGenPt[i]         = new TH2D(Form("hEleMu_RecoGenPt%s", Sample[i]),"",6,2.5,14.5,6,2.5,14.5);
    hEleMu_RecoGenTauPt[i]         = new TH2D(Form("hEleMu_RecoGenTauPt%s", Sample[i]),"",6,2.5,14.5,12,2.5,12.5);
    hEleMu_RecoGenInvmass[i]          = new TH2D(Form("hEleMu_RecoGenInvmass%s",Sample[i]),"",5,4,24,5,4,24);
    hEleMu_RecoGenvSumPt[i]       = new TH2D(Form("hEleMu_RecoGenvSumPt%s",Sample[i]),"",5,0,10,5,0,10);




    cout << "file " << TauTau[i]->GetEntries()  << endl;
    ReadTauTauTree  TauTauR(TauTau[i]);
    TauTauR.fChain->SetBranchStatus("*",1);
    if (TauTauR.fChain == 0) return;
    Long64_t nentries = TauTauR.fChain->GetEntriesFast();
    cout << TauTau[i]->GetName() << "    " << nentries << endl;
    
    Long64_t nbytes = 0, nb = 0;
    for (Long64_t jentry=0; jentry<nentries;jentry++) {
   // for (Long64_t jentry=0; jentry<100;jentry++) {
      Long64_t ientry_evt_R = TauTauR.LoadTree(jentry);
      if (ientry_evt_R < 0) break;
      nb = TauTauR.fChain->GetEntry(jentry);   nbytes += nb;
     

      // "===================      Fill gen Info here ======================     
      //if(TauTauR.EleMuEvents == 0 ) continue;
      if(TauTauR.hasMuon == 0 ) continue;
      if(TauTauR.Check == 0 ) continue;
       if(TauTauR.DeltaR_mu > 0.15) continue;
      if(abs(TauTauR.Tau_Eta2) > 2.4 ) continue;

     double ET = TauTauR.Gen_Mu_Eta2;
     double PT = TauTauR.Gen_Mu_Pt2;
  //   double Weight1 = tnp_weight_softid_upc_pbpb( PT, ET, 0);
    // double Weight2 = tnp_weight_trigger_upc_pbpb( PT, ET, 0);
//     double W = Weight1*Weight2*lumiNormGUPC;
    // double W = Weight1*Weight2*lumiNormSC;
//     double W = 1;

    
      hTau_GenPt[i]->Fill(TauTauR.Tau_Pt2,W);
 
      // ===================   Apply analysis cuts on reco MC ======================           
      bool isPassing =// TauTauR.Ele_Pt > 2.0 && abs(TauTauR.Ele_Eta) <  2.2
 			 abs(TauTauR.Mu_Eta) < 2.4 && TauTauR.Mu_Pt > 2.5;	
      double ET2 = TauTauR.Mu_Eta;
      double PT2 = TauTauR.Mu_Pt;
      //double Weight3 = tnp_weight_softid_upc_pbpb( PT2, ET2, 0);
      //double Weight4 = tnp_weight_trigger_upc_pbpb( PT2, ET2, 0);
      //double W2 = Weight3*Weight4*lumiNormGUPC;
     // double W2 = Weight3*Weight4*lumiNormSC;
//      double W2 = 1;


			
      if(TauTauR.ok_trigger == 1) { //trigger 
	if(isPassing){
//        cout << "Weight3:" << Weight3 << endl;
       	
	hMu_Pt[i]->Fill(TauTauR.Mu_Pt,W2);
/* 
        hEleMu_RecoGenTauPt[i]->Fill(TauTauR.Mu_Pt,TauTauR.Tau_Pt2,W2);
        hEleMu_RecoGenPt[i]->Fill(TauTauR.Mu_Pt,TauTauR.Gen_Mu_Pt2,W2);
        hEleMu_RecoGenInvmass[i]->Fill(TauTauR.vSum_M,TauTauR.Gen_vSum_M2,W2);
        hEleMu_RecoGenvSumPt[i]->Fill(TauTauR.vSum_Pt,TauTauR.Gen_vSum_Pt2,W2);
*/	
        hEleMu_RecoGenTauPt[i]->Fill(TauTauR.Tau_Pt2,TauTauR.Mu_Pt,W2);


	} //isPassing analysis cuts
      } //trigger
      
      
    } //entry
 // } // for 4 files
  
/*************************************************************/
TFile *f2=new TFile("RootFile/GenMatch/SC_25thSep_mass10bin.root","r");
//TFile *f1=new TFile("RootFile/GenMatch/GAMMA_UPC_Hist.root","r");
//TH1D *hMu_Pt_SC = (TH1D*)f2->Get("muon_pt_all_TAUTAU");
//TH1D *hEleMu_Mass_SC = (TH1D*)f2->Get("dilepton_mass_all_TAUTAU");
//TH1D *hEleMu_vSumPt_SC = (TH1D*)f2->Get("dilepton_pt_all_TAUTAU");
/*************************************************************/




  /***************************  Start reading data and subtract background contributions ********************************/

// TFile *f1=new TFile("RootFile/EleMu_histograms_withSFs_withTag.root","r");
// TFile *f1=new TFile("RootFile/GenMatch/Data_MuMuFSR_19thSep.root","r");
 TFile *f1=new TFile("RootFile/GenMatch/Data_MuMuFSR_4thOctober.root","r");
//  TFile *f1=new TFile("RootFile/Data_MuMuFSR_19thSep_Mass20bin.root");
 TH1D *hMu_Pt_data    = (TH1D*)f1->Get("muon_pt_all_Data");
 TH1D *hDilepton_Mass_data = (TH1D*)f1->Get("dilepton_mass_all_Data");
 TH1D *hDilepton_Pt_data = (TH1D*)f1->Get("dilepton_pt_all_Data");

 TH1D *hMu_Pt_mumufsr    = (TH1D*)f1->Get("muon_pt_all_MUMU_FSR");
 TH1D *hDilepton_Mass_mumufsr = (TH1D*)f1->Get("dilepton_mass_all_MUMU_FSR");
 TH1D *hDilepton_Pt_mumufsr = (TH1D*)f1->Get("dilepton_pt_all_MUMU_FSR");

 TH1D* hMu_Pt_data_bkgSub       = subtractBkg(hMu_Pt_data, hMu_Pt_mumufsr, "hMu_Pt_data_bkgSub", 6,2.5,14.5);
 TH1D* hDilepton_Mass_data_bkgSub       = subtractBkg(hDilepton_Mass_data, hDilepton_Mass_mumufsr, "hDilepton_Mass_data_bkgSub", 5,4,24);
 TH1D* hDilepton_Pt_data_bkgSub       = subtractBkg(hDilepton_Pt_data, hDilepton_Pt_mumufsr, "hDilepton_Pt_data_bkgSub", 5,0,10);

// "===================     Unfold data ======================     
/*********************For Non-regularization*******************************/
/*
      TUnfold::ERegMode regMode=TUnfold::kRegModeNone;
      TUnfold::EConstraint constraintMode=TUnfold::kEConstraintNone;
      TUnfoldDensity::EDensityMode densityFlags = TUnfoldDensity::kDensityModeBinWidth;
      TUnfold::EHistMap outputMap = TUnfold::kHistMapOutputHoriz;

   TUnfoldDensity unfold_TauPt(hEleMu_RecoGenTauPt[0],outputMap,regMode,constraintMode,densityFlags);
 //  TUnfoldDensity unfold_MuPt(hEleMu_RecoGenPt[0],outputMap,regMode,constraintMode,densityFlags);
   //TUnfoldDensity unfold_Mass(hEleMu_RecoGenInvmass[0],outputMap,regMode,constraintMode,densityFlags);
   //TUnfoldDensity unfold_vSumPt(hEleMu_RecoGenvSumPt[0],outputMap,regMode,constraintMode,densityFlags);

   double tau = -50.0;

   unfold_TauPt.DoUnfold(tau, hMu_Pt_data_bkgSub);
  // unfold_MuPt.DoUnfold(tau, hMu_Pt_data_bkgSub);
  // unfold_Mass.DoUnfold(tau, hDilepton_Mass_data_bkgSub);
  // unfold_vSumPt.DoUnfold(tau,hDilepton_Pt_data_bkgSub);

   TH1 *hUnfold_TauPt_data=unfold_TauPt.GetOutput("");
  // TH1 *hUnfold_MuPt_data=unfold_MuPt.GetOutput("");
  // TH1 *hUnfold_Mass_data=unfold_Mass.GetOutput("");
  // TH1 *hUnfold_VSumPt_data=unfold_vSumPt.GetOutput("");
*/
/********************************For regularization***************************/

  TUnfold::ERegMode regMode=TUnfold::kRegModeSize;


  TUnfold::EConstraint constraintMode=TUnfold::kEConstraintNone;
  TUnfoldDensity::EDensityMode densityFlags = TUnfoldDensity::kDensityModeBinWidth;
  TUnfold::EHistMap outputMap = TUnfold::kHistMapOutputHoriz;
  TUnfold unfold(hEleMu_RecoGenTauPt[0],
  TUnfold::kHistMapOutputHoriz,
  TUnfold::kRegModeSize,
  TUnfold::kEConstraintNone);
  unfold.SetInput( hMu_Pt_data_bkgSub,1.0);
  TGraph *lcurve=0;
  TSpline *logTauX=0,*logTauY=0;
  cout << "Unfold1" << endl; 
  unfold.ScanLcurve(20,0.,0.,&lcurve,&logTauX,&logTauY);
  cout << "Unfold2" << endl;
  double tau=unfold.GetTau();
  double logTau=TMath::Log10(tau);
  double lcurveX=logTauX->Eval(logTau);
  double lcurveY=logTauY->Eval(logTau);


  unfold.DoUnfold(tau,hMu_Pt_data_bkgSub);
  TH1 *hUnfold_TauPt_data;
  unfold.GetOutput(hUnfold_TauPt_data);


/************************** ********************************/
  
double GenNorm;
double RecoNorm;

if(GenScale){
   GenNorm = hMu_Pt_data_bkgSub->Integral(0,6)/hTau_GenPt[0]->Integral(0,6);
   hTau_GenPt[0]->Scale(GenNorm);
}

if(RecoScale){
   RecoNorm = hMu_Pt_data_bkgSub->Integral(0,6)/hMu_Pt[0]->Integral(0,6);
   hMu_Pt[0]->Scale(RecoNorm);
}


/*
 // cout << " invariant mass reco " << Sample[0] <<  " :" << hEleMu_Mass[0]->Integral() << endl;  
//  cout << " invariant mass Data" <<   " :" << hDilepton_Mass_data->Integral() << endl;  
  cout << " invariant mass mumuFSR" <<   " :" << hDilepton_Mass_mumufsr->Integral() << endl;  
  cout << " invariant mass Data substract bkg" <<   " :" << hDilepton_Mass_data_bkgSub->Integral() << endl;  
  cout << " invariant mass Gen " <<   " :" << hEleMu_Gen_Mass[0]->Integral() << endl;  
  cout << " Pt Gen Mu" <<   " :" << hMu_Gen_Pt[0]->Integral() << endl;  
//  cout << " Pt Gen Ele" <<   " :" << hEle_Gen_Pt[0]->Integral() << endl;  
  cout << " Pt Gen EleMu" <<   " :" << hEleMu_Gen_vSumPt[0]->Integral() << endl;  
  cout << "Pt Data bkgSub:" << ":"  << hMu_Pt_data_bkgSub->Integral() << endl;
  cout << "VSumPt:" << ":" << hDilepton_Pt_data_bkgSub->Integral() << endl;
*/
  outf->cd();
  
  //outf->Write();

  hUnfold_TauPt_data->Write();

  hMu_Pt_data_bkgSub->Write();
  /* 
  hUnfoldTauPt->Write();
  hUnfoldPt->Write();
  hUnfoldMass->Write();
  hUnfoldVSumPt->Write();
*/

  hMu_Pt[0]->Write();

  hTau_GenPt[0]->Write();

/****************************************Draw**********************************/
  /************************************TauPt*****************************/
//  TCanvas*c8 =  new TCanvas();
  TCanvas* c8 = new TCanvas("c8","TauPt",50,50,600,670);
  TPad *pad7 = new TPad("pad7", "pad7", 0, 0.3, 1, 1.0);
  pad7->SetBottomMargin(0);
  pad7->Draw();
  pad7->cd();
  hUnfold_TauPt_data->GetXaxis()->SetTitle("p_{T} of Muon (Green #tau)");
  hUnfold_TauPt_data->GetYaxis()->SetTitle("Events");
  hUnfold_TauPt_data->SetLineColor(kBlack);
  hUnfold_TauPt_data->SetMarkerStyle(21);
  hUnfold_TauPt_data->GetYaxis()->SetRangeUser(-15, 40);
  hUnfold_TauPt_data->Draw("p");
  
  hMu_Pt[0]->SetLineColor(kBlue);
  hMu_Pt[0]->GetXaxis()->SetTitle("p_{T} of Muon (Green #tau)");
  hMu_Pt[0]->GetYaxis()->SetTitle("Events");
 // hMu_Pt[0]->SetMarkerStyle(22);
 // hMu_Pt[0]->SetMarkerColor(9);
  hMu_Pt[0]->SetLineWidth(2);
  hMu_Pt[0]->Draw("histsame");


/*  hUnfoldTauPt->SetLineColor(kRed);
  hUnfoldTauPt->GetXaxis()->SetTitle("p_{T} of Muon (Green #tau)");
  hUnfoldTauPt->GetYaxis()->SetTitle("Events");
  hUnfoldTauPt->Draw("p");
  hUnfold_TauPt_data->GetXaxis()->SetTitle("p_{T} of Muon (Green #tau)");
  hUnfold_TauPt_data->GetYaxis()->SetTitle("Events");
  hUnfold_TauPt_data->SetLineColor(kBlack);
  hUnfold_TauPt_data->SetMarkerStyle(21);
  hUnfold_TauPt_data->Draw("psame");*/

  hMu_Pt[0]->SetLineColor(kBlue);
  //hMu_Pt[0]->Draw("histsame");
  hTau_GenPt[0]->SetLineColor(kGreen);
  //hTau_GenPt[0]->SetMarkerStyle(23);
  //hTau_GenPt[0]->SetMarkerColor(3);
  hTau_GenPt[0]->SetLineWidth(2);
  hTau_GenPt[0]->Draw("histsame");
  hMu_Pt_data_bkgSub->SetLineColor(kViolet);
 // hMu_Pt_data_bkgSub->SetMarkerStyle(20);
 // hMu_Pt_data_bkgSub->SetMarkerColor(6);
  hMu_Pt_data_bkgSub->SetLineWidth(2);
  hMu_Pt_data_bkgSub->Draw("histsame");
  TLegend *leg4 = new TLegend(0.60, 0.72, 0.80, 0.9);
  leg4->SetFillColor(0);
  leg4->SetBorderSize(0);
  leg4->SetFillStyle(0);
  leg4->SetTextFont(43);
  leg4->SetTextSize(15);
//  leg4->AddEntry(hUnfoldTauPt,"Unfolded MC","pl");
  leg4->AddEntry(hTau_GenPt[0], "Gen-level MC","pl");
  leg4->AddEntry(hMu_Pt[0], "Reconstructed MC", "pl");
  leg4->AddEntry(hMu_Pt_data_bkgSub, "Data", "pl");
  leg4->AddEntry(hUnfold_TauPt_data, "Unfolded Data", "pl");
  leg4->Draw();
  c8->cd();
  TPad *pad8 = new TPad("pad8", "pad8", 0, 0.06, 1, 0.28);
  pad8->SetTopMargin(0);
  pad8->SetBottomMargin(0.3);
  pad8->SetGridx();
  pad8->Draw();
  pad8->cd();
  //c6->cd();
  TString strRatio4=hTau_GenPt[0]->GetName();
  strRatio4+="_over_";
  strRatio4+=hTau_GenPt[0]->GetName();
  TH1D* hratio4=(TH1D*)hUnfold_TauPt_data->Clone(strRatio4);
  hratio4->Sumw2();
  hratio4->Divide(hTau_GenPt[0]);
  hratio4->SetTitle("");
  hratio4->GetYaxis()->SetTitle("UnfoData/Gen");
  hratio4->SetLabelSize(0.1, "XYZ");
  hratio4->SetLabelFont(42, "XYZ");
  hratio4->SetLabelOffset(0.007, "XYZ");
  hratio4->SetTitleSize(0.11, "XYZ");
  hratio4->GetXaxis()->SetTitleOffset(1.0);
  hratio4->GetYaxis()->SetTitleOffset(0.3);
  hratio4->SetMarkerStyle(21);
  hratio4->GetYaxis()->SetRangeUser(0, 2);
  hratio4->Draw("p");
  c8->SaveAs("Plot/Data/TauPt_Data.png");
/****************************************************************/
   outf->Close();

  
}

TH1D* subtractBkg(TH1D *hdata, TH1D *hmumufsr, const char* name, int nbins, float xmin, float xmax){

  TH1D *hans = new TH1D(name,"",nbins,xmin,xmax);
  //TH1D* hans = (TH1D*)hdata->Clone();
   cout << " Name of histogram :" << name << endl;
   for (int i=1; i<=hdata->GetNbinsX(); i++) {
      double binCont  = hdata->GetBinContent(i)-hmumufsr->GetBinContent(i);
      double binError = sqrt(pow(hdata->GetBinError(i),2) + pow(hmumufsr->GetBinError(i),2) );


      hans->SetBinContent(i, binCont);
      hans->SetBinError(i, binError);
      //cout << "******bin:" << i << " " << hdata->GetBinContent(i) << " " << hqed->GetBinContent(i) << " " << hcep->GetBinContent(i) << "   " <<  hans->GetBinContent(i) << "+/-" <<  binError<< endl;
      //double err=0;
      //for (int ivar=1; ivar<14; ivar++) err += pow(hvari[ivar]->GetBinContent(i)-hvari[0]->GetBinContent(i),2);
      //hans->SetBinError(i,sqrt(err+pow(hans->GetBinError(i),2)+pow(glob_syst,2)));
   }

   return hans;
}




TH1D* getXSecHist(TH1D *hist, const char* name, int nbins, float xmin, float xmax, float luminosity){

  TH1D *hans = new TH1D(name,"",nbins,xmin,xmax);
   cout << " Name of histogram :" << name << endl;
   for (int i=1; i<=hist->GetNbinsX(); i++) {

      double binCont  = hist->GetBinContent(i)/(luminosity*hist->GetBinWidth(i));

      double binError = binCont * sqrt(pow(hist->GetBinError(i)/hist->GetBinContent(i),2));
      //double binError = hist->GetBinError(i);
      if(binCont == 0){ hans->SetBinContent(i, 0);
       hans->SetBinError(i, 0);}
	else{
      hans->SetBinContent(i, binCont);
      hans->SetBinError(i, binError);
    }

      //cout << "******bin:" << i << " " << hans->GetBinContent(i) << "+/-" <<  binError<< endl;
      //double err=0;
      //for (int ivar=1; ivar<14; ivar++) err += pow(hvari[ivar]->GetBinContent(i)-hvari[0]->GetBinContent(i),2);
      //hans->SetBinError(i,sqrt(err+pow(hans->GetBinError(i),2)+pow(glob_syst,2)));
   }
   return hans;

}


void printOutput(TH1D *hist){
 for (int i=1; i<=hist->GetNbinsX(); i++) {
      cout << "******bin:" << i << " " << hist->GetBinContent(i) << "+/-" <<  hist->GetBinError(i) << endl; 
    }
}
void make_canvas(TCanvas *& canvas){

  int W = 600;
  int H = 670;
  float T = 0.08;
  float B = 0.14; 
  float L = 0.18;
  float R = 0.04;

  canvas->SetFillColor(0);
  canvas->SetBorderMode(0);
  canvas->SetFrameFillStyle(0);
  canvas->SetFrameBorderMode(0);
  canvas->SetLeftMargin( L );
  canvas->SetRightMargin( R );
  canvas->SetTopMargin( T );
  canvas->SetBottomMargin( B );
  canvas->SetTickx(0);
  canvas->SetTicky(0);
}

void make_canvas_ratio(TCanvas *& can){

  int W = 700;
  int H = 600;
  float T = 0.08;
  float B = 0.14; 
  float L = 0.14;
  float R = 0.04;

  can->SetFillColor(0);
  can->SetBorderMode(0);
  can->SetFrameFillStyle(0);
  can->SetFrameBorderMode(0);
  can->SetLeftMargin( L );
  can->SetRightMargin( R );
  can->SetTopMargin( T );
  can->SetBottomMargin( B );
  can->SetTickx(0);
  can->SetTicky(0);
}


void make_hist(TH1D *& hist, Color_t kcolor, int kstyle){
  // For the axis labels:
  hist->SetLabelColor(1, "XYZ");
  hist->SetLabelFont(42, "XYZ");
  hist->SetLabelOffset(0.007, "XYZ");
  hist->SetLabelSize(0.05, "XYZ");
 
  hist->SetMarkerColor(kcolor);
  hist->SetFillColor(kcolor);
  hist->SetMarkerStyle(kstyle);
 
}

void make_hist_ratio(TH1D *& hist, Color_t kcolor, int kstyle){
  // For the axis labels:
  hist->SetLabelColor(1, "XYZ");
  hist->SetLabelFont(42, "XYZ");
  hist->SetLabelOffset(0.007, "XYZ");
  hist->SetLabelSize(0.1, "XYZ");

  // For the axis titles:
  hist->SetTitleFont(42, "XYZ");
  hist->SetTitleSize(0.11, "XYZ");
  hist->GetXaxis()->SetTitleOffset(1.1);
  hist->GetYaxis()->SetTitleOffset(0.4);
  hist->SetMarkerColor(kcolor);
  hist->SetFillColor(kcolor);
  hist->SetMarkerStyle(kstyle);
 
}


float getAcoBinSF(float acop) {
      if (acop<0.002) return 0.990761;
      else if (acop<0.004) return 0.962675;
      else if (acop<0.006) return 0.861166;
      else if (acop<0.008) return 0.598683;
      else if (acop<0.01) return 0.264046;    
      else return 1;
   
}