C++ TH1F::Divide方法代码示例

//--- function to divide histograms then draw---
void drawDivHist(const char* hn1, const char* hn2, const char* opt, const char* name, const char* title, const int nbin, const float min, const float max, bool log=false, const int lc=0, const int ls=0, const int lw=0)
{
   // find input histos
   if (!TString(opt).Contains("same")) printf("\n");
   printf("h1: %s. h2: %s. Draw: %s\n", hn1, hn2, name);
   TH1F * h1; TH1F * h2;
   if (gROOT->FindObject(hn1))
      h1 = dynamic_cast<TH1F*>(gROOT->FindObject(hn1));
   else {
      printf("%s is not found, please check the histogram name\n", hn1);
      return;
   }
   if (gROOT->FindObject(hn2))
      h2 = dynamic_cast<TH1F*>(gROOT->FindObject(hn2));
   else {
      printf("%s is not found, please check the histogram name\n", hn2);
      return;
   }

   //--- Make/set histogram ---
   printf("hist: %s %d %f %f\n",name,nbin,min,max);
   TH1F * h = createHist(name, title, nbin, min, max);
   setHist(h,lc,ls,lw);

   //--- Action ---
   h->Divide(h1,h2);

   //--- Draw ---
   TCanvas * c = makeCanvas(name,title,log,opt);
   h->Draw(opt);
}

//--------------------------------------------------------------------------------------------------
// make fake rate histogram from numerator and denominators
//--------------------------------------------------------------------------------------------------
TH1F* createFakeRatePtHist(vector<string> datasetFiles, vector<string> datasetNames,
                           int faketype, string histName) {

  string numeratorHistName = "";
  string denominatorHistName = "";
  string dirName = "";

  if (faketype == 11) {
    dirName = "ComputeElectronFakeRateMod";
    numeratorHistName = "hElectronNumeratorEt";
    denominatorHistName = "hElectronDenominatorEt";
  } else if (faketype == 13) {
    dirName = "ComputeMuonFakeRateMod";
    numeratorHistName = "hMuonNumeratorPt";
    denominatorHistName = "hMuonDenominatorPt";
  } else {
    cerr << "Error: faketype = " << faketype << " is not recognized." << endl;
    return 0;
  }

  TH1F *denominator = addAllSamples(datasetFiles, datasetNames, dirName, denominatorHistName);
  TH1F *numerator = addAllSamples(datasetFiles, datasetNames, dirName, numeratorHistName);

  //create fake rate hist from denominator binning
  TH1F *fakeRateHist = (TH1F*)numerator->Clone(histName.c_str());
  fakeRateHist->GetYaxis()->SetTitle("Fake Rate");
  //divide by the denominator to get fake rate;
  fakeRateHist->Divide(denominator);

  return fakeRateHist;
}

void make_histos_syst_rawyield(TString file_syst, TString file_default, TString out_tag){

  TFile *f1 = new TFile(file_syst.Data(),"read");
  TFile *f2 = new TFile(file_default.Data(),"read");
  TDirectoryFile *dir1 = (TDirectoryFile*)(f1->Get("effunf"));
  TDirectoryFile *dir2 = (TDirectoryFile*)(f2->Get("effunf"));
  TList *list = dir1->GetListOfKeys();

  TFile *f = new TFile(Form("plots/ratiosyst_%s.root",out_tag.Data()),"recreate");

  for (int i=0; i<list->GetSize(); i++){
    TString name = dir1->GetListOfKeys()->At(i)->GetName();
    if (!(name.Contains("hreco_"))) continue;
    TObject *obj1 = dir1->Get(name.Data());
    assert(obj1);
    TObject *obj2 = dir2->Get(name.Data());
    assert(obj2);
    TString newname = name;
    newname.Append("_ratiosyst");
    if (name.EndsWith("_0")) newname.ReplaceAll("_0_","_EBEB_");
    if (name.EndsWith("_1")) newname.ReplaceAll("_1_","_EBEE_");
    if (name.EndsWith("_2")) newname.ReplaceAll("_2_","_EEEE_");
    TH1F *h = (TH1F*)(((TH1F*)obj1)->Clone(newname.Data()));
    h->SetTitle(h->GetName());
    h->Divide((TH1F*)obj2);
    for (int j=0; j<h->GetNbinsX(); j++) h->SetBinError(j+1,0);
    for (int j=0; j<h->GetNbinsX(); j++) h->SetBinContent(j+1,1+fabs(1-h->GetBinContent(j+1)));
    f->cd();
    h->Write();
  }

}

static void	make_plot( const char* i )
{
	c.cd(++pad_id);

	char	*hist_name = Form("h%s_hit", i);
	TH1F	*hit = new TH1F(hist_name, "", time_bins, min_time, max_time);

	events->Draw(
		Form("timestamp >> %s", hist_name),
		Form("(t%s_hits_count == 4) && (t%s_track_count >= 1)", i, i)
		);

	hist_name = Form("h%s_all", i);
	TH1F	*all = new TH1F(hist_name, "", time_bins, min_time, max_time);

	events->Draw(
		Form("timestamp >> %s", hist_name)
		);

	hist_name = Form("h%s_eff", i);
	TH1F	*eff = new TH1F(hist_name, "", time_bins, min_time, max_time);
	eff->Divide(hit, all);
	eff->GetYaxis()->SetRangeUser(0.0, 1.0);
	eff->Draw();

	c.cd(++pad_id);
	hist_name = Form("h%s_eff_distrib", i);
	TH1F	*eff_distrib = new TH1F(hist_name, "", 100, 0, 0.5);
	int nbins = eff->GetXaxis()->GetNbins();
	for (int bin = 1; bin <= nbins; bin++) {
		eff_distrib->Fill(eff->GetBinContent(bin));
	}
	eff_distrib->Draw();
}

void FindZ2Weight(){

  NCanvas(1,1,"data");
  NCanvas(1,1,"ratio");

  TH1F * Pyt;
  TH1F * Dat;
  TH1F * Data;

  TFile *fzee = new TFile("ZDiffOutputfile.root");
  Pyt  =  (TH1F*)fzee->Get("NoCuts_GoodVtx_Z2PY6");
  Dat  =  (TH1F*)fzee->Get("NoCuts_GoodVtx_DATA10");
  Data  =  (TH1F*)fzee->Get("NoCuts_GoodVtx_DATA10");
  data->cd(1);
  Data->Draw();

  Float_t scale = Pyt->Integral()/Dat->Integral();
  cout << "Scaling factor to conserve MC luminosity = " << scale << endl;
  ratio->cd(1);
  Dat->Divide(Pyt);
  //  Dat->Scale(scale);
  Dat->Draw();
  cout << "Float_t Z2Weight[] = {" ;
  for (Int_t j = 1;j<=19; j++)   cout <<  Dat->GetBinContent(j) << ", " ;
  cout << " 0 }; " << endl;
}

void beffAnalysis(const char* input, const char* output) {
  // input
  TFile* inputFile = TFile::Open(input);
  TTree* btagEff = (TTree*)inputFile->Get("btagEff");

  // output
  TFile* outputFile = TFile::Open(output,"RECREATE");

  // histogram with proper binning... cloned later on 
  Double_t binning[23] = {20,25,30,35,40,45,50,55,60,65,70,75,80,85,90,95,100,120,140,160,180,200,1000};
  TH1F* ptSpectrum = new TH1F("PtSpectrum","PtSpectrum",22,binning);
  ptSpectrum->Sumw2();

  // produce the ratio plot for the 12 combinations of (CSVL,CSVM,CSVT),(Barrel,Endcap),(b,c,l)
  TClonesArray algorithms("TCut",3);
  new(algorithms[0]) TCut("CSVL","csv>0.244");
  new(algorithms[1]) TCut("CSVM","csv>0.679");
  new(algorithms[2]) TCut("CSVT","csv>0.898");

  TClonesArray etaRegions("TCut",2);
  new(etaRegions[0]) TCut("Barrel","abs(eta)<=1.2");
  new(etaRegions[1]) TCut("Endcaps","abs(eta)>1.2");

  TClonesArray flavor("TCut",3);
  new(flavor[0]) TCut("l","abs(flavor)!=4 && abs(flavor)!=5");
  new(flavor[1]) TCut("c","abs(flavor)==4");
  new(flavor[2]) TCut("b","abs(flavor)==5");

  for(int i=0; i< algorithms.GetEntries() ; ++i) {
    outputFile->mkdir(((TCut*)algorithms.At(i))->GetName());
    outputFile->cd(((TCut*)algorithms.At(i))->GetName());
    for(int j=0; j< etaRegions.GetEntries() ; ++j) {
      for(int k=0; k< flavor.GetEntries() ; ++k) {
        // histogram before tagging
        TH1F* pretag = ptSpectrum->Clone("pretag");
        btagEff->Draw("pt>>pretag",((*(TCut*)etaRegions.At(j))&&(*(TCut*)flavor.At(k)))*"eventWeight");
        // histogram after tagging
        TH1F* posttag = ptSpectrum->Clone("posttag");
        btagEff->Draw("pt>>posttag",((*(TCut*)algorithms.At(i))&&(*(TCut*)etaRegions.At(j))&&(*(TCut*)flavor.At(k)))*"eventWeight");
        // ratio
        TH1F* ratio = posttag->Clone(Form("h_eff_bTagOverGoodJet_pt%s_%s",((TCut*)flavor.At(k))->GetName(),
                                                                          ((TCut*)etaRegions.At(j))->GetName()));
        ratio->Divide(pretag);
        // cleanup
        delete pretag;
        delete posttag;
      }
    }
  }

  // cleanup
  algorithms.Delete();
  etaRegions.Delete();
  flavor.Delete();
  ptSpectrum->SetDirectory(0);
  outputFile->Write();
  outputFile->Close();
  inputFile->Close();
}

void reweight_eta_1d(TH1F* weight_eta,TH1F* weight_etao,TH2F*weight_etan,TH2F* weight_eta2o,TH2F* weight_etanr,TH2F*weight_eta2,RooDataSet **dset, RooDataSet *dsetdestination, int numvar){

  if (!(*dset)) return;

  TH1F *hnum = new TH1F("hnum","hnum",n_etabins_forreweighting,etabins_forreweighting);
 TH1F *hden = new TH1F("hden","hden",n_etabins_forreweighting,etabins_forreweighting);
//  TH1F *hnum = new TH1F("hnum","hnum",25,0.,2.5);
//  TH1F *hden = new TH1F("hden","hden",25,0.,2.5);
  hnum->Sumw2();
  hden->Sumw2();

  const char* etaname=Form("rooeta%d",numvar);

  for (int i=0; i<(*dset)->numEntries(); i++){
    hden->Fill(fabs((*dset)->get(i)->getRealValue(etaname)),(*dset)->store()->weight(i));
  }
  for (int i=0; i<dsetdestination->numEntries(); i++){
    hnum->Fill(fabs(dsetdestination->get(i)->getRealValue(etaname)),dsetdestination->store()->weight(i));
  }
  hnum->Scale(1.0/hnum->Integral());
  hden->Scale(1.0/hden->Integral());

  hnum->Divide(hden);
  TH1F *h = hnum;

  RooDataSet *newdset = new RooDataSet(**dset,Form("%s_etarew",(*dset)->GetName()));
  newdset->reset();
  for (int i=0; i<(*dset)->numEntries(); i++){
    RooArgSet args = *((*dset)->get(i));
    float oldw = (*dset)->store()->weight(i);
    float eta = args.getRealValue(etaname);
    float neww = oldw*h->GetBinContent(h->FindBin(fabs(eta)));
    if(debug){ 
		weight_eta->Fill(neww);	
		weight_etao->Fill(oldw);
		weight_etan->Fill(h->FindBin(fabs(eta)),neww);	
		weight_eta2o->Fill(h->FindBin(fabs(eta)),oldw);
		if(oldw!=0 && neww!=0)weight_etanr->Fill(h->FindBin(fabs(eta)),oldw/neww);
		else {weight_etanr->Fill(-10,1);}
	   // weight_pt2->Fill(pt,neww/oldw);
	   if(oldw!=0 && neww!=0)weight_eta2->Fill(fabs(eta),oldw/neww);
	   else {weight_eta2->Fill(-10,1);}
	}  
	  newdset->add(args,neww);
	  
	  }


  newdset->SetName((*dset)->GetName());
  newdset->SetTitle((*dset)->GetTitle());
  delete hnum; delete hden;

  RooDataSet *old_dset = *dset;
  *dset=newdset;
  std::cout << "Eta 1d rew: norm from " << old_dset->sumEntries() << " to " << newdset->sumEntries() << std::endl;

  delete old_dset;

};

void PlotMinEtFromSim(Bool_t isPhos = kFALSE){
  gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);
  Float_t min = 0;
  float max = 1;
  TString filename, detname;
  if(isPhos){
    min = 0.655;
    max = 0.785;
    detname = "PHOS";
    filename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.PHOS.LHC11a10a_bis.Run139465.root";
  }
  else{
    min = 0.58;
    max = 0.725;
    filename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal.LHC11a10a_bis.Run139465.root";
    detname = "EMCal";
  }
  
  TFile *f = TFile::Open(filename, "READ");
  TList *l = dynamic_cast<TList*>(f->Get("out1"));
  TH1F *fHistSimulatedGammaEnergyAboveThreshold = l->FindObject("fHistSimulatedGammaEnergyAboveThreshold");
  TH1F *fHistSimulatedGammaEnergy = l->FindObject("fHistSimulatedGammaEnergy");
  SetStyles(fHistSimulatedGammaEnergyAboveThreshold,20,TColor::kRed);
  fHistSimulatedGammaEnergyAboveThreshold->Divide(fHistSimulatedGammaEnergy);

    TCanvas *c1 = new TCanvas("c1","Simulation",600,400);
    c1->SetTopMargin(0.02);
    c1->SetRightMargin(0.03);
    c1->SetLeftMargin(0.11745);
    c1->SetBottomMargin(0.11745);
    c1->SetBorderSize(0);
    c1->SetFillColor(0);
    c1->SetFillColor(0);
    c1->SetBorderMode(0);
    c1->SetFrameFillColor(0);
    c1->SetFrameBorderMode(0);
    fHistSimulatedGammaEnergyAboveThreshold->SetMaximum(max +0.1);
    fHistSimulatedGammaEnergyAboveThreshold->SetMinimum(min-0.1);
    fHistSimulatedGammaEnergyAboveThreshold->GetXaxis()->SetTitle("Centrality bin");
    fHistSimulatedGammaEnergyAboveThreshold->GetYaxis()->SetTitle("f_{minEt}");
    fHistSimulatedGammaEnergyAboveThreshold->GetYaxis()->SetLabelSize(0.06);
    fHistSimulatedGammaEnergyAboveThreshold->GetXaxis()->SetLabelSize(0.06);
    fHistSimulatedGammaEnergyAboveThreshold->GetYaxis()->SetTitleSize(0.06);
    fHistSimulatedGammaEnergyAboveThreshold->GetXaxis()->SetTitleSize(0.06);
    fHistSimulatedGammaEnergyAboveThreshold->Draw();
    TLine *lineMin = new TLine(-0.5,min,19.5,min);
    lineMin->Draw();
    TLine *lineMax = new TLine(-0.5,max,19.5,max);
    lineMax->Draw();
    lineMin->SetLineColor(TColor::kBlue);
    lineMax->SetLineColor(TColor::kBlue);
    lineMin->SetLineStyle(2);
    lineMax->SetLineStyle(2);

    TString outfile = "/tmp/MinEtFromSim"+detname+".png";
    c1->SaveAs(outfile.Data());
}

//--------------------------------------------------------------------------------------------------
// make fake rate histogram from numerator and denominators
//--------------------------------------------------------------------------------------------------
// TGraphAsymmErrors* createFakeRatePtHist(vector<string> datasetFiles, vector<string> datasetNames, 
//                            int denominatortype, int faketype, int chargetype,  
//                            int eventType, string histName) {
TH1F* createFakeRatePtHist(vector<string> datasetFiles, vector<string> datasetNames, 
                           int denominatortype, int faketype, int chargetype,  
                           int eventType, string histName) {

  string numeratorHistName = "";
  string denominatorHistName = "";
  string dirName = "";

  string chargeTypeName = "";
  if (chargetype == -1)
    chargeTypeName = "Minus";
  else if (chargetype == 1)
    chargeTypeName = "Plus";

  string eventTypeName = "";
  if (eventType == -1)
    eventTypeName = "MinusW";
  else if (eventType == 1)
    eventTypeName = "PlusW";

  string denominatorTypeName = "";
  if (denominatortype == 1)
    denominatorTypeName = "Track";

  if (faketype == 11) {
    dirName = "ComputeElectronFakeRateMod";
    numeratorHistName = "h" + chargeTypeName + "Electron" + eventTypeName + "NumeratorPt";
    denominatorHistName = "h" + chargeTypeName + "Electron" + eventTypeName + denominatorTypeName + "DenominatorPt";
  } else if (faketype == 13) {
    dirName = "ComputeMuonFakeRateMod";
    numeratorHistName = "h" + chargeTypeName + "Muon" + eventTypeName + "NumeratorPt";
    denominatorHistName = "h" + chargeTypeName + "Muon" + eventTypeName + denominatorTypeName + "DenominatorPt";
  } else {
    cerr << "Error in createFakeRatePtHist: faketype = " << faketype << " is not recognized." << endl;
    assert(false);
  }

  TH1F *denominator = addAllSamples(datasetFiles, datasetNames, dirName, denominatorHistName);
  TH1F *numerator = addAllSamples(datasetFiles, datasetNames, dirName, numeratorHistName);

  //create fake rate hist from denominator binning
  TH1F *fakeRateHist = (TH1F*)numerator->Clone(histName.c_str());
  fakeRateHist->GetYaxis()->SetTitle("Fake Rate");
  //divide by the denominator to get fake rate;
  fakeRateHist->Divide(denominator);

  //TGraphAsymmErrors* fakeRate = (numerator, denominator);
  

  //for (int i=1 ; i <= fakeRateHist->GetXaxis()->GetNbins(); i++) {
    //cout << i << " : " << fakeRateHist->GetXaxis()->GetBinCenter(i) << " " << fakeRateHist->GetBinContent(i) << "  +- " << fakeRateHist->GetBinError(i) << endl;
    //if (fakeRateHist->GetBinContent(i) > 0 && fakeRateHist->GetBinError(i) < 0.01) {
      //fakeRateHist->SetBinError(i,0.01);
    //}
  //}

  return fakeRateHist;
}

void reweight_diphotonpt_1d(TH1F*weight_diphopt,TH1F*weight_diphopto,TH2F*weight_diphoptn,TH2F*weight_diphopt2o,TH2F*weight_diphoptnr, TH2F* weight_diphopt2,RooDataSet **dset, RooDataSet *dsetdestination){

  if (!(*dset)) return;
  if (!(dsetdestination)) return;
///numerator and denominator
  TH1F *hnum = new TH1F("hnum","hnum",n_diphoptbins_forreweighting,diphoptbins_forreweighting);
  TH1F *hden = new TH1F("hden","hden",n_diphoptbins_forreweighting,diphoptbins_forreweighting);
  hnum->Sumw2();
  hden->Sumw2();

     
// RooAbsData->get*() Load a given row of data

 //getRealValue Get value of a RooAbsReal stored in set with given name. If none is found, value of defVal is returned.
 for (int i=0; i<(*dset)->numEntries(); i++){
      hden->Fill(fabs((*dset)->get(i)->getRealValue("roodiphopt")),(*dset)->store()->weight(i));
  }
  for (int i=0; i<dsetdestination->numEntries(); i++){
    hnum->Fill(fabs(dsetdestination->get(i)->getRealValue("roodiphopt")),dsetdestination->store()->weight(i));
  }
 //normalize to one 
  hnum->Scale(1.0/hnum->Integral());
  hden->Scale(1.0/hden->Integral());
  hnum->Divide(hden);
  TH1F *h = hnum;
  RooDataSet *newdset = new RooDataSet(**dset,Form("%s_diphoptrew",(*dset)->GetName()));
  newdset->reset();
  assert(newdset);
  for (int i=0; i<(*dset)->numEntries(); i++){
    RooArgSet args = *((*dset)->get(i));
    float oldw = (*dset)->store()->weight(i);
    float diphopt = fabs((*dset)->get(i)->getRealValue("roodiphopt"));
    float neww = oldw*h->GetBinContent(h->FindBin(fabs(diphopt)));
    newdset->add(args,neww);
    if(debug){
		weight_diphopt->Fill(neww);
		weight_diphopto->Fill(oldw);
		weight_diphopt2o->Fill((h->FindBin(fabs(diphopt))),oldw);	
		weight_diphoptn->Fill((h->FindBin(fabs(diphopt))),neww);	
		if(oldw!=0)weight_diphoptnr->Fill(h->FindBin(fabs(diphopt)),oldw/neww);
		else {weight_diphoptnr->Fill(-10,1);}//cout << "dipho weight old 0" << endl;}
		if(oldw!=0)weight_diphopt2->Fill(diphopt,oldw/neww);
		else {weight_diphopt2->Fill(-10,1);}//cout << "dipho weight old 0" << endl;}
	 }
		}

  newdset->SetName((*dset)->GetName());
  newdset->SetTitle((*dset)->GetTitle());
  delete hnum; delete hden;
 
  RooDataSet *old_dset = *dset;
  *dset=newdset;
  std::cout << "Diphoton Pt 1d rew: norm from " << old_dset->sumEntries() << " to " << newdset->sumEntries() << std::endl;

  delete old_dset;

};

void draw_ratio(std::vector<TH1F*> h,
	  TString name, TString xTitle,
	  double xmin, double xmax,
	  TString legHeader = "", bool legRIGHT = true, bool legTOP = true,
	  bool logX = false, bool stat = false, int rebin = -1, int orbin = -1,
	  TString option = "", int nclus = 99) {  //double ymin_ratio, double ymax_ratio,
  TCanvas* can = new TCanvas(name+"_ratio",name+"_ratio",900,450);
  can->cd();

  double legxmin = (legRIGHT ? 0.55 : 0.18);
  double legxmax = legxmin+0.25;
  double legymin = (legTOP ? 0.70 : 0.15);
  double legymax = legymin+0.15;
  TLegend* leg = new TLegend(legxmin,legymin,legxmax,legymax);
  if (legHeader!="") leg->SetHeader(legHeader);
  leg->SetTextSize(0.04);
  leg->SetFillColor(0);
  leg->SetLineColor(0);

  TString options = (option=="" ? "pe" : option);

  if (rebin>0) h[h.size()-1]->Rebin(rebin); //to rebin benchmark before divide
  for (size_t i=0; i<h.size(); i++) {
    //if(h[i]->GetNbinsX() != orbin) cout << "WARNING: orbin for " << h[i]->GetName() << " are " << h[i]->GetNbinsX() << endl; //debug - shift of h[][] wrt clu[][]
    if (rebin>0 && i<(h.size()-1)) h[i]->Rebin(rebin);
    TH1F* ratio = (TH1F*)h[i]->Clone("ratio_"+name);
    ratio->Sumw2();
    ratio->Divide(h[h.size()-1]); //benchmark is at the end.

    if (logX) gPad->SetLogx();
    ratio->SetMarkerStyle(20+i);
    ratio->SetMarkerSize(1.0); //1.2
    ratio->GetXaxis()->SetRangeUser(xmin,xmax);
    ratio->SetMinimum(-0.1);
    ratio->SetMaximum(4);
    if (i==0){    //just for the first one
     ratio->GetXaxis()->SetLabelSize(0.05);
     ratio->GetXaxis()->SetTitle(xTitle);
     ratio->GetXaxis()->SetTitleOffset(1);
     ratio->GetXaxis()->SetTitleSize(0.06);
     ratio->GetYaxis()->SetTitle("ratio");
     ratio->GetYaxis()->SetTitleSize(0.06);
     ratio->GetYaxis()->SetTitleOffset(0.7);
     ratio->GetYaxis()->SetLabelSize(0.05);
    }
    string nam = "";
    nam = translate(clu[nclus][0].c_str());
    if(i==(h.size()-1)) leg->AddEntry(h[i],nam.c_str(),"l"); //to print only benchmark (first in the list)
    if (i==1) options = options + (stat ? "sames" : "same"); //once is enought
    ratio->Draw(options);
  }
  leg->Draw("same");
  drawPrivate(0.04);
  can->Update();
  can->SaveAs(Outfolder+name+"_ratio.png");
}

void findZscale(){
  TFile* file1 = new TFile("DoubleEG_prompt.root");
  TFile* file2 = new TFile("nob_pt_zelel_cnt_bkg.root");
  TH1F* hist1 = (TH1F*)file1->Get("ana/nob_pt_zelel_cnt");
  TH1F* hist2 = (TH1F*)file2->Get("total bkg");
  TH1F* divided = (TH1F*)hist1->Clone();
  divided->Divide(hist1, hist2);
  gStyle->SetOptStat(0000);
  divided->GetYaxis()->SetTitle("Data/MC");
  divided->Draw("ep");
  divided->SaveAs("ElectronFit.root");
}

void compare_Zsub_4253(){

  TFile *f42 = new TFile("histo_xsec_dphi_inclusive_DPHI_SCALE_07_42.root");
  TH1F *h42 = (TH1F*)(f42->Get("xsec"));
  
  TFile *f42up = new TFile("histo_xsec_dphi_inclusive_DPHI_SCALE_1_PURITY_1_42.root");
  TH1F *h42up = (TH1F*)(f42up->Get("xsec"));

  TFile *f42unf = new TFile("histo_finalxs_fortheorycomp_dphi_42.root");
  TH1F *h42unf = (TH1F*)(f42unf->Get("histo_finalxs_fortheorycomp_dphi"));
  h42unf->SetMarkerStyle(1);
  h42unf->SetLineColor(kBlack);
  h42unf->SetFillColor(kBlack);
  h42unf->SetMarkerColor(kBlack);
  TH1F *h42unforig = (TH1F*)(h42unf->Clone("histo_finalxs_fortheorycomp_dphi_orig"));

  TFile *f53unf = new TFile("histo_finalxs_fortheorycomp_dphi.root");
  TH1F *h53unf = (TH1F*)(f53unf->Get("histo_finalxs_fortheorycomp_dphi"));
  h53unf->SetMarkerStyle(1);
  h53unf->SetLineColor(kRed);
  h53unf->SetFillColor(kRed);
  h53unf->SetMarkerColor(kRed);
  TH1F *h53unforig = (TH1F*)(h53unf->Clone("histo_finalxs_fortheorycomp_dphi_orig"));

  TFile *f53unc = new TFile("histo_systsummaryfinal_individual_dphi_inclusive.root");
  TH1F *h53unc = (TH1F*)(f53unc->Get("systplot_zee_allcat"));
  
  TH1F *ratio = (TH1F*)(h42->Clone("ratio42"));
  ratio->Divide(h42up);

  for (int bin=0; bin<h42->GetNbinsX(); bin++){
    float x = h42unf->GetBinContent(bin+1);
    float r = ratio->GetBinContent(bin+1);
    h42unf->SetBinError(bin+1,x*fabs(r-1));
    cout << r << endl;
  }

  for (int bin=0; bin<h53unf->GetNbinsX(); bin++){
    float x = h53unf->GetBinContent(bin+1);
    float e = h53unc->GetBinContent(bin+1);
    h53unf->SetBinError(bin+1,x*e);
  }

  h42unforig->SetFillStyle(3004);
  h53unforig->SetFillStyle(3005);

  h42unforig->Draw("E2");
  h53unforig->Draw("E2 same");
  h42unf->Draw("E1 same");
  h53unf->Draw("E1 same");


}

void ZVXT1HF0()
{

  NCanvas(1,1,"data");
  NCanvas(1,1,"ratio");

  TH1F * HDatHF0;
  TH1F * HDat;

  TFile *fzee = new TFile("ZDiffOutputfile_DataOnly_2010_2011.root");
  // TFile *fzee = new TFile("ZDiffOutputfile.root");
  HDat  =  (TH1F*)fzee->Get("NVTX1_InstLumiPerBx_DATA10");
  
  HDatHF0  =  (TH1F*)fzee->Get("HF0NVTX1_InstLumiPerBx_DATA10");
  data->cd(1);
  HDat->Draw();
  NSetTitle(HDat,"Luminosity [10^{-30} cm^{-2} s^{-1}]", "Entries");
  HDat->SetTitle("Number of Z produced in events without PU");
  HDatHF0->Draw("SAME HIST");
  NHSetMarker(HDat,2,20,0.4);
  TLegend *legend = new TLegend(0.6,0.7,0.9,0.8);
  legend->SetTextFont(72);
  legend->SetTextSize(0.03);
  legend->SetBorderSize(0);


  legend->AddEntry(HDat,"All Z","p");
  legend->AddEntry(HDatHF0,"Z with minE_HF=0","l");
  legend->Draw();


  TH1F * Hra = (TH1F * ) HDatHF0->Clone();
  //  Hra->GetXaxis->SetRange(0,12);
  NStat(Hra,0);
  Hra->Divide(HDat);
  NLogAxis(0,1);
  ratio->cd(1);

  Hra->Draw();
  Hra->SetTitle("Fraction of Z with (PU=0 && E_HF=0)");
  Hra->Fit("expo","","",0.,1.2);
  NSetTitle(Hra,"Luminosity [10^{-30} cm^{-2} s^{-1}]", "Fraction");
  NLogAxis(0,1);
  NText(.05,0.1, "0.62");
  NText(0.3,0.1, "2.47");
  NText(0.7,0.1, "4.93");
  NText(1.1,0.1, "7.4");
  NText(1.4,0.1, "<Number of Int.>");


}

void Ef_Branch(){


  gROOT->ProcessLine(" .L tdrstyle.C");
  setTDRStyle();
  gStyle->SetOptStat(0);

  TH1F *IsoPaftl = new TH1F("IsoPaftl","IsoP cut only",100,0,100); 
  TH1F *IsoPaftlb = new TH1F("IsoPaftlb","IsoP cut only",100,0,100); 
  

  //Pt dependency

  IsoPaftl->Divide(EffNp,EffN0,1.,1.,"B");
  IsoPaftlb->Divide(EffNpb,EffN0b,1.,1.,"B");


  TCanvas *c10 = new TCanvas("c10","Branch Out Cuts",900,600);
  c10->Divide(3,2);

  IsoPaftl->Draw();


  IsoPaftlb->Draw("same");

  IsoPaftl->SetLineColor(kRed);
  IsoPaftl->SetMarkerColor(kRed);


  IsoPaftl->SetMarkerSize(0.5);
  IsoPaftl->GetYaxis()->SetTitle("Only iso P Cut Efficiency");
  IsoPaftl->GetXaxis()->SetTitle("Nvtx");
  
  c10->Update();
  c10->SaveAs("LCheck.png");


}