#include <iostream>
#include <string>
void reproduceAN(){
  TFile *myFile = TFile::Open("/home/christos/public_html/Pythia/pythia8226/examples/hist1.root");
  TTree *tree = (TTree*)myFile->Get("tree");

  string acceptanceCut = "((Reco_QQ_mumi_4mom.Pt()>3.5&&TMath::Abs(Reco_QQ_mumi_4mom.Eta())<1.2)||(Reco_QQ_mumi_4mom.Pt()>(5.77-1.8*TMath::Abs(Reco_QQ_mumi_4mom.Eta()))&&TMath::Abs(Reco_QQ_mumi_4mom.Eta())>=1.2&&TMath::Abs(Reco_QQ_mumi_4mom.Eta())<2.1)||(Reco_QQ_mumi_4mom.Pt()>1.8&&TMath::Abs(Reco_QQ_mumi_4mom.Eta())>=2.1&&TMath::Abs(Reco_QQ_mumi_4mom.Eta())<2.4))&&((Reco_QQ_mupl_4mom.Pt()>3.5&&TMath::Abs(Reco_QQ_mupl_4mom.Eta())<1.2)||(Reco_QQ_mupl_4mom.Pt()>(5.77-1.8*TMath::Abs(Reco_QQ_mupl_4mom.Eta()))&&TMath::Abs(Reco_QQ_mupl_4mom.Eta())>=1.2&&TMath::Abs(Reco_QQ_mupl_4mom.Eta())<2.1)||(Reco_QQ_mupl_4mom.Pt()>1.8&&TMath::Abs(Reco_QQ_mupl_4mom.Eta())>=2.1&&TMath::Abs(Reco_QQ_mupl_4mom.Eta())<2.4))"; 
  string analysisCut = "Reco_QQ_mupl_nPixWMea>0&&Reco_QQ_mupl_nTrkWMea>5&&Reco_QQ_mupl_dxy<.3&&Reco_QQ_mupl_dz<20&&Reco_QQ_mupl_dz>-20&&Reco_QQ_mupl_isGoodMuon&&Reco_QQ_mumi_nPixWMea>0&&Reco_QQ_mumi_nTrkWMea>5&&Reco_QQ_mumi_dxy<.3&&Reco_QQ_mumi_dz<20&&Reco_QQ_mumi_dz>-20&&Reco_QQ_mumi_isGoodMuon&&Reco_QQ_VtxProb>.01";
  string ppDecayCut = "1";
  string massChargeCut = "Reco_QQ_4mom.M()>2.946&&Reco_QQ_4mom.M()<3.346&&Reco_QQ_sign==0";
  string matchCut = "TMath::Sqrt(TMath::Power(Reco_QQ_mupl_4mom.Eta()-Gen_QQ_mupl_4mom.Eta(),2)+TMath::Power(Reco_QQ_mupl_4mom.Phi()-Gen_QQ_mupl_4mom.Phi(),2))<.03&&TMath::Sqrt(TMath::Power(Reco_QQ_mumi_4mom.Eta()-Gen_QQ_mumi_4mom.Eta(),2)+TMath::Power(Reco_QQ_mumi_4mom.Phi()-Gen_QQ_mumi_4mom.Phi(),2))<.03";
  string otherCut_num = "TMath::Abs(Reco_mu_4mom.Rapidity())<1.6";
  string otherCut_denom = "TMath::Abs(Gen_mu_4mom.Rapidity())<1.6";
  
  TH1D *histoPt_denom = new TH1D("histoPt_denom","P_{T}",4,0,12);
  histoPt_denom->GetXaxis()->SetTitle("P_{T} (GeV)");

  TH1D *histoPt_num = new TH1D("histoPt_num","P_{T} (After Cuts)",4,0,12);
  histoPt_num->GetXaxis()->SetTitle("P_{T} (GeV)");

  TH1D *histoPt_eff = new TH1D("histoPt_eff","P_{T} Eff.",4,0,12);
  histoPt_eff->GetXaxis()->SetTitle("P_{T} (GeV)");
  histoPt_eff->GetYaxis()->SetRange(0,1);
  
  const char *denom_cut = (acceptanceCut+"&&"+otherCut_denom).c_str();
  const char *num_cut = (ppDecayCut+"&&"+massChargeCut+"&&"+matchCut+"&&"+otherCut_num+"&&"+analysisCut).c_str();

  TCanvas *c0 = new TCanvas();
  myTree.Draw("Gen_mu_4mom.Pt()>>histoPt_denom",denom_cut,"COLZ");

  TCanvas *c1 = new TCanvas();
  myTree.Draw("Reco_mu_4mom.Pt()>>histoPt_num",num_cut,"COLZ");

  histoPt_eff->Sumw2();
  histoPt_eff->Divide(histoPt_num,histoPt_denom,1,1,"B");

  TCanvas *c2 = new TCanvas();
  histoPt_eff->SetMaximum(1);
  histoPt_eff->SetMinimum(0);
  histoPt_eff->Draw("COLZ");
}