C++ TH1D::SetMarkerStyle方法代码示例

//.........这里部分代码省略.........
    SeeddeltaEtaforSimTrackvalidHist->GetXaxis()->SetTitle("(recEta-simEta)");
    SeeddeltaEtaforSimTrackvalidHist->GetXaxis()->CenterTitle();
    SeeddeltaEtaforSimTrackvalidHist->Draw();
    string SeeddeltaEtaforSimTrackvalidCanvasName = FinalOutput + "SeeddeltaEtaforSimTrackvalid" + OutputPlotNameFix;
    SeeddeltaEtaforSimTrackvalidCanvas->SaveAs(SeeddeltaEtaforSimTrackvalidCanvasName.c_str());


    TCanvas* SeedPurityforSimTrackvalidCanvas = new TCanvas("SeedPurityforSimTrackvalidCanvas", "SeedPurityforSimTrackvalidCanvas", 800, 600);
    SeedPurityforSimTrackvalidCanvas->cd();
    SeedPurityforSimTrackvalidHist->Draw();
    string SeedPurityforSimTrackvalidCanvasName = FinalOutput + "SeedPurityforSimTrackvalid" + OutputPlotNameFix;
    SeedPurityforSimTrackvalidCanvas->SaveAs(SeedPurityforSimTrackvalidCanvasName.c_str());

    TCanvas* ChargeCheckforSimTrackvalidCanvas = new TCanvas("ChargeCheckforSimTrackvalidCanvas", "ChargeCheckforSimTrackvalidCanvas", 800, 600);
    ChargeCheckforSimTrackvalidCanvas->cd();
    double HistEntries = ChargeCheckforSimTrackvalidHist->GetEntries() / 100.;
    ChargeCheckforSimTrackvalidHist->Scale(1./HistEntries);
    ChargeCheckforSimTrackvalidHist->GetXaxis()->SetTitle("simCharge*recCharge");
    ChargeCheckforSimTrackvalidHist->GetXaxis()->CenterTitle(1);
    ChargeCheckforSimTrackvalidHist->GetYaxis()->SetTitle("fraction %");
    ChargeCheckforSimTrackvalidHist->GetYaxis()->CenterTitle(1);
    ChargeCheckforSimTrackvalidHist->Draw();
    string ChargeCheckforSimTrackvalidCanvasName = FinalOutput + "ChargeCheckforSimTrackvalid" + OutputPlotNameFix;
    ChargeCheckforSimTrackvalidCanvas->SaveAs(ChargeCheckforSimTrackvalidCanvasName.c_str());

    TCanvas* SeedNumberforSimTrackvalidCanvas = new TCanvas("SeedNumberforSimTrackvalidCanvas", "SeedNumberforSimTrackvalidCanvas", 800, 600);
    SeedNumberforSimTrackvalidCanvas->cd();
    SeedNumberforSimTrackvalidHist->Draw();
    string SeedNumberforSimTrackvalidCanvasName = FinalOutput + "SeedNumberforSimTrackvalid" + OutputPlotNameFix;
    SeedNumberforSimTrackvalidCanvas->SaveAs(SeedNumberforSimTrackvalidCanvasName.c_str());

    TCanvas* SeedEfficiencyforSimTrackvalidCanvas = new TCanvas("SeedEfficiencyforSimTrackvalidCanvas", "SeedEfficiencyforSimTrackvalidCanvas", 800, 600);
    SeedEfficiencyforSimTrackvalidCanvas->cd();
    SeedEfficiencyforSimTrackvalidHist->Draw();
    string SeedEfficiencyforSimTrackvalidCanvasName = FinalOutput + "SeedEfficiencyforSimTrackvalid" + OutputPlotNameFix;
    SeedEfficiencyforSimTrackvalidCanvas->SaveAs(SeedEfficiencyforSimTrackvalidCanvasName.c_str());

    TCanvas* SeedNumberforSimTrackinvalidCanvas = new TCanvas("SeedNumberforSimTrackinvalidCanvas", "SeedNumberforSimTrackinvalidCanvas", 800, 600);
    SeedNumberforSimTrackinvalidCanvas->cd();
    SeedNumberforSimTrackinvalidHist->Draw();
    string SeedNumberforSimTrackinvalidCanvasName = FinalOutput + "SeedNumberforSimTrackinvalid" + OutputPlotNameFix;
    SeedNumberforSimTrackinvalidCanvas->SaveAs(SeedNumberforSimTrackinvalidCanvasName.c_str());

    TCanvas* SeedEfficiencyforSimTrackinvalidCanvas = new TCanvas("SeedEfficiencyforSimTrackinvalidCanvas", "SeedEfficiencyforSimTrackinvalidCanvas", 800, 600);
    SeedEfficiencyforSimTrackinvalidCanvas->cd();
    SeedEfficiencyforSimTrackinvalidHist->Draw();
    string SeedEfficiencyforSimTrackinvalidCanvasName = FinalOutput + "SeedEfficiencyforSimTrackinvalid" + OutputPlotNameFix;
    SeedEfficiencyforSimTrackinvalidCanvas->SaveAs(SeedEfficiencyforSimTrackinvalidCanvasName.c_str());

    double SeedEfficiencyforSimTrackinvalid = 100. * SeedEfficiencyforSimTrackinvalidHist->GetMean();
    double SeedEfficiencyforSimTrackvalid = 100. * SeedEfficiencyforSimTrackvalidHist->GetMean();
    SeedEfficiencyHist->SetBinContent(1, SeedEfficiencyforSimTrackinvalid);
    SeedEfficiencyHist->SetBinContent(2, SeedEfficiencyforSimTrackvalid);
    SeedEfficiencyHist->GetXaxis()->SetBinLabel(1, "for invalid simTrack");
    SeedEfficiencyHist->GetXaxis()->SetBinLabel(2, "for valid simTrack");
    TCanvas* SeedEfficiencyCanvas = new TCanvas("SeedEfficiencyCanvas", "SeedEfficiencyCanvas", 800, 600);
    SeedEfficiencyCanvas->cd();
    SeedEfficiencyHist->GetYaxis()->SetTitle("Efficiency %");
    SeedEfficiencyHist->GetYaxis()->CenterTitle(1);
    SeedEfficiencyHist->SetMarkerStyle(3);
    SeedEfficiencyHist->SetMarkerSize(3);
    SeedEfficiencyHist->Draw("P");
    string SeedEfficiencyCanvasName = FinalOutput + "SeedEfficiency" + OutputPlotNameFix;
    SeedEfficiencyCanvas->SaveAs(SeedEfficiencyCanvasName.c_str());

    TCanvas* RecBendingPhi2PtCanvas = new TCanvas("RecBendingPhi2PtCanvas", "RecBendingPhi2PtCanvas", 800, 600);
    RecBendingPhi2PtCanvas->cd();
    RecBendingPhi2PtHist->Draw();
    string RecBendingPhi2PtCanvasName = FinalOutput + "RecBendingPhi2Pt" + OutputPlotNameFix;
    RecBendingPhi2PtCanvas->SaveAs(RecBendingPhi2PtCanvasName.c_str());

    TCanvas* PtRatoofRecBendingPhiCanvas = new TCanvas("PtRatoofRecBendingPhiCanvas", "PtRatoofRecBendingPhiCanvas", 800, 600);
    PtRatoofRecBendingPhiCanvas->cd();
    PtRatoofRecBendingPhiHist->Draw();
    string PtRatoofRecBendingPhiCanvasName = FinalOutput + "PtRatoofRecBendingPhi" + OutputPlotNameFix;
    PtRatoofRecBendingPhiCanvas->SaveAs(PtRatoofRecBendingPhiCanvasName.c_str());

    TCanvas* RecBendingLastPhiCanvas = new TCanvas("RecBendingLastPhiCanvas", "RecBendingLastPhiCanvas", 800, 600);
    RecBendingLastPhiCanvas->cd();
    RecBendingLastPhiHist->Draw();
    string RecBendingLastPhiCanvasName = FinalOutput + "RecBendingLastPhi" + OutputPlotNameFix;
    RecBendingLastPhiCanvas->SaveAs(RecBendingLastPhiCanvasName.c_str());

    Int_t linsav = gStyle->GetLineWidth();
    gStyle->SetLineWidth(2);
    TCanvas* SimReverseBendingCanvas = new TCanvas("SimReverseBendingCanvas", "SimReverseBendingCanvas", 800, 800);
    SimReverseBendingCanvas->cd();
    TPad* SimReverseBendingPad = new TPad("SimReverseBendingPad", "SimReverseBendingPad", 0, 0, 1, 1);
    SimReverseBendingPad->Draw();
    SimReverseBendingPad->cd();
    SimReverseBendingPad->Range(-800, -800, 800, 800);
    unsigned int segmentNumber = SimReverseBending->GetEntries();
    cout << "Number of segments: " << segmentNumber << endl;
    for(unsigned int j = 0; j < segmentNumber; j++) {
        ((TLine*)(SimReverseBending->At(j)))->Print();
        ((TLine*)(SimReverseBending->At(j)))->Draw("SAME");
    }
    string SimReverseBendingCanvasName = FinalOutput + "SimReverseBending" + OutputPlotNameFix;
    SimReverseBendingCanvas->SaveAs(SimReverseBendingCanvasName.c_str());
}

TF1 *fit(TTree *nt,TTree *ntMC,double ptmin,double ptmax, bool ispPb, int count){   
   //cout<<cut.Data()<<endl;
   //static int count=0;
   //count++;
   TCanvas *c= new TCanvas(Form("c%d",count),"",600,600);
   TH1D *h = new TH1D(Form("h%d",count),"",50,5,6);
   TH1D *hMC = new TH1D(Form("hMC%d",count),"",50,5,6);

   TString iNP="7.26667e+00*Gaus(x,5.10472e+00,2.63158e-02)/(sqrt(2*3.14159)*2.63158e-02)+4.99089e+01*Gaus(x,4.96473e+00,9.56645e-02)/(sqrt(2*3.14159)*9.56645e-02)+3.94417e-01*(3.74282e+01*Gaus(x,5.34796e+00,3.11510e-02)+1.14713e+01*Gaus(x,5.42190e+00,1.00544e-01))";
   TF1 *f = new TF1(Form("f%d",count),"[0]*([7]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+(1-[7])*Gaus(x,[1],[8])/(sqrt(2*3.14159)*[8]))+[3]+[4]*x+[5]*("+iNP+")");
   nt->Project(Form("h%d",count),"mass",Form("%s&&pt>%f&&pt<%f",seldata_2y.Data(),ptmin,ptmax));   
   ntMC->Project(Form("hMC%d",count),"mass",Form("%s&&pt>%f&&pt<%f",seldata_2y.Data(),ptmin,ptmax));   
   clean0(h);

   TH1D *hraw = new TH1D(Form("hraw%d",count),"",50,5,6);
   clean0(hraw);
   hraw = (TH1D*)h->Clone(Form("hraw%d",count));


   h->Draw();
   f->SetParLimits(4,-1000,0);
   f->SetParLimits(2,0.01,0.05);
   f->SetParLimits(8,0.01,0.05);
   f->SetParLimits(7,0,1);
   f->SetParLimits(5,0,1000);

   f->SetParameter(0,setparam0);
   f->SetParameter(1,setparam1);
   f->SetParameter(2,setparam2);
   f->SetParameter(8,setparam3);
   f->FixParameter(1,fixparam1);
   h->GetEntries();

   hMC->Fit(Form("f%d",count),"q","",5,6);
   hMC->Fit(Form("f%d",count),"q","",5,6);
   f->ReleaseParameter(1);
   hMC->Fit(Form("f%d",count),"L q","",5,6);
   hMC->Fit(Form("f%d",count),"L q","",5,6);
   hMC->Fit(Form("f%d",count),"L q","",5,6);
   hMC->Fit(Form("f%d",count),"L m","",5,6);

   f->FixParameter(1,f->GetParameter(1));
   f->FixParameter(2,f->GetParameter(2));
   f->FixParameter(7,f->GetParameter(7));
   f->FixParameter(8,f->GetParameter(8));
   
   h->Fit(Form("f%d",count),"q","",5,6);
   h->Fit(Form("f%d",count),"q","",5,6);
   f->ReleaseParameter(1);
   h->Fit(Form("f%d",count),"L q","",5,6);
   h->Fit(Form("f%d",count),"L q","",5,6);
   h->Fit(Form("f%d",count),"L q","",5,6);
   h->Fit(Form("f%d",count),"L m","",5,6);
   h->SetMarkerSize(0.8);
   h->SetMarkerStyle(20);
   cout <<h->GetEntries()<<endl;

   // function for background shape plotting. take the fit result from f
   TF1 *background = new TF1(Form("background%d",count),"[0]+[1]*x");
   background->SetParameter(0,f->GetParameter(3));
   background->SetParameter(1,f->GetParameter(4));
   background->SetParameter(2,f->GetParameter(5));
   background->SetParameter(3,f->GetParameter(6));
   background->SetLineColor(4);
   background->SetRange(5,6);
   background->SetLineStyle(2);
   
   // function for signal shape plotting. take the fit result from f
   TF1 *Bkpi = new TF1(Form("fBkpi%d",count),"[0]*("+iNP+")");
   Bkpi->SetParameter(0,f->GetParameter(5));
   Bkpi->SetLineColor(kGreen+1);
   Bkpi->SetFillColor(kGreen+1);
//   Bkpi->SetRange(5.00,5.28);
   Bkpi->SetRange(5.00,6.00);
   Bkpi->SetLineStyle(1);
   Bkpi->SetFillStyle(3004);

   // function for signal shape plotting. take the fit result from f
   TF1 *mass = new TF1(Form("fmass%d",count),"[0]*([3]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+(1-[3])*Gaus(x,[1],[4])/(sqrt(2*3.14159)*[4]))");
   mass->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2),f->GetParameter(7),f->GetParameter(8));
   mass->SetParError(0,f->GetParError(0));
   mass->SetParError(1,f->GetParError(1));
   mass->SetParError(2,f->GetParError(2));
   mass->SetParError(7,f->GetParError(7));
   mass->SetParError(8,f->GetParError(8));
   mass->SetLineColor(2);
   mass->SetLineStyle(2);

//   cout <<mass->Integral(0,1.2)<<" "<<mass->IntegralError(0,1.2)<<endl;
   h->SetMarkerStyle(24);
   h->SetStats(0);
   h->Draw("e");
   h->SetXTitle("M_{B} (GeV/c^{2})");
   h->SetYTitle("Entries / (20 MeV/c^{2})");
   h->GetXaxis()->CenterTitle();
   h->GetYaxis()->CenterTitle();
   h->SetTitleOffset(1.5,"Y");
   h->SetAxisRange(0,h->GetMaximum()*1.2,"Y");
   Bkpi->Draw("same");
   background->Draw("same");   
//.........这里部分代码省略.........

void Zselection() {
  TLatex *tplus = labelLatex(0.20,0.80,"7 TeV Data");
  gROOT->SetStyle("Plain"); 
  gStyle->SetOptStat(0);
  gStyle->SetPalette(1);
  //the correct W x-sec to use is 31314, c.f. 24170, 30380
  double w_scale = 35.0;//(31314.0 / 24170.0) * 3.2;
  double z_scale = w_scale;
  double qcd_scale = 14./1000.0;
  double data_scale = 1.;
  unsigned int rbin = 5;
  unsigned int numPlots = 11;
  bool doPrint = false;

  TString folder = "hltmu15_goodevsel/Zselection";

  TFile *file0 = TFile::Open("results/" + folder + "/RecoRoutines_Z-selection_realdata.root");
  TFile *file1 = TFile::Open("results/" + folder + "/RecoRoutines_Z-selection_WJets_madgraph_June2010.root");
  //TFile *file1 = TFile::Open("results/" + folder + "/RecoRoutines_W-selection_WJets_sherpa.root");
  //TFile *file2 = TFile::Open("results/" + folder + "/RecoRoutines_W-selection_QCD_AllPtBins_7TeV_Pythia.root");
  TFile *file3 = TFile::Open("results/" + folder + "/RecoRoutines_Z-selection_ZJets_madgraph_June2010.root");
  //TFile *file4 = TFile::Open("results/MuPt10/RecoRoutines_W-selection_TTbarJets_tauola_madgraph_June2010.root");

  //dynamic array size not allowed in CINT hehe
  TCanvas *canvas [11] = {makeCanvas("reco_wpt_plus"), makeCanvas("reco_wpt_minus"), makeCanvas("reco_wpt"), makeCanvas("muon_pt_plus"), makeCanvas("muon_pt_minus"), makeCanvas("muon_pt"), makeCanvas("pf_mt_plus"), makeCanvas("pf_mt_minus"), makeCanvas("pf_mt"), makeCanvas("lpvar_plus"), makeCanvas("lpvar_minus")};
  TLegend *legend [11] = {makeLegend(), makeLegend(), makeLegend(), makeLegend(), makeLegend(), makeLegend(), makeLegend(), makeLegend(), makeLegend(), makeLegend(), makeLegend()};

  TString plotnames [11] = {"RECO_PolPlots_50toinf/RECO_pfMHTPlus",
                            "RECO_PolPlots_50toinf/RECO_pfMHTMinus",
                            "RECO_PolPlots_50toinf/RECO_pfMHT", 
                            "RECO_PolPlots_50toinf/RECO_MuonPtPlus", 
                            "RECO_PolPlots_50toinf/RECO_MuonPtMinus", 
                            "RECO_PolPlots_50toinf/RECO_MuonPt", 
                            "RECO_PolPlots_50toinf/RECO_pfMTPlus",
                            "RECO_PolPlots_50toinf/RECO_pfMTMinus",
                            "RECO_PolPlots_50toinf/RECO_pfMT", 
                            "RECO_PolPlots_50toinf/RECO_ICVarPlus", 
                            "RECO_PolPlots_50toinf/RECO_ICVarMinus"};
  TString plotlabels [11] = {"RECO P_{T}(W+) [GeV]", "RECO P_{T}(W-) [GeV]", "RECO P_{T}(W) [GeV]", "P_{T}(Muon+) [GeV]", "P_{T}(Muon-) [GeV]", "P_{T}(Muon) [GeV]", "M_{T}+ [GeV]", "M_{T}- [GeV]", "M_{T} [GeV]", "LP(+)", "LP(-)"};
  int rbinextra [11] = {2,2,2,4,4,4,2,2,2,4,4};
  double xmin [11] = {50., 50., 50., 0., 0., 0., 30., 30., 30., -0.5, -0.5};
  double xmax [11] = {120., 120., 120., 140., 140., 140., 90., 90., 90., 1.5, 1.5};
  double ymax [11] = {170., 150., 300., 150., 100., 250., 160., 130., 260., 150., 80.};
  int skip [11] = {0,0,0,0,0,0,0,0,0,0,0};


  for(unsigned int i=0; i<numPlots; i++) {
    if(skip[i] == 0) {
    TH1D * data = (TH1D*)file0->Get(plotnames[i]);
    TH1D * w = (TH1D*)file1->Get(plotnames[i]);
    //TH1D * qcd = (TH1D*)file2->Get(plotnames[i]);
    TH1D * z = (TH1D*)file3->Get(plotnames[i]);
    //TH1D * tt = (TH1D*)file4->Get(plotnames[i]);

    w->Rebin(rbin*rbinextra[i]);
    //qcd->Rebin(rbin*rbinextra[i]);
    z->Rebin(rbin*rbinextra[i]);
    //tt->Rebin(rbin);
    data->Scale(data_scale);
    data->Rebin(rbin*rbinextra[i]);

    w->Scale(w_scale);
    z->Scale(z_scale);
    //tt->Scale(lumi_scale);
    //qcd->Scale(qcd_scale);

    TH1D * mc = (TH1D*)w->Clone();
    //mc->Add(w);
    mc->Add(z);
    //mc->Add(tt);
    mc->SetLineWidth(3);
    //mc->SetLineStyle(9);
    mc->SetLineColor(kGray);
    //qcd->SetLineColor(kGreen);
    //qcd->SetLineWidth(2);
    //qcd->GetXaxis()->SetRangeUser(xmin[i],xmax[i]);
    //qcd->GetXaxis()->SetTitle(plotlabels[i]);
    //qcd->GetXaxis()->SetTitleSize(0.05);
    //qcd->GetYaxis()->SetRangeUser(0.1,ymax[i]);
    //qcd->GetYaxis()->SetTitle("Events / 2.8 pb^{-1}");
    //qcd->GetYaxis()->SetTitleSize(0.05);
    //qcd->GetYaxis()->SetTitleOffset(1.36);
    canvas[i]->cd(1);//->SetLogy();
    //qcd->DrawCopy("h");
    w->SetLineColor(kBlue);
    w->SetLineWidth(2);
    w->DrawCopy("h");
    z->SetLineColor(kRed);
    z->SetLineWidth(2);
    z->DrawCopy("sameh");
    //tt->SetLineColor(kGold);
    //tt->SetLineWidth(2);
    //tt->DrawCopy("sameh");

    mc->DrawCopy("sameh");
    data->SetMarkerStyle(20);
    data->DrawCopy("samep");
    //tplus->DrawClone("same");
    legend[i]->AddEntry(data, "7TeV Data", "p");
    legend[i]->AddEntry(mc, "All MC", "l");
//.........这里部分代码省略.........

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// drawFigures4and5
//
// Zpt
// LeadingJetPt
// Njets
//
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
void drawFigures4and5(TString var = "Zpt") 
{
  gSystem->mkdir("pdf", kTRUE);
  gSystem->mkdir("png", kTRUE);

  gInterpreter->ExecuteMacro("WZPaperStyle.C");

  Bool_t drawMcfm = (var.EqualTo("Zpt")) ? true : false;

  TString variable;
  TString xtitle;
  TString ytitle;

  if (var.EqualTo("Zpt"))          variable = "Z";
  if (var.EqualTo("LeadingJetPt")) variable = "leading jet";
  if (var.EqualTo("Njets"))        variable = "N_{jets}";

  if (var.EqualTo("Njets"))
    {
      xtitle = variable;
      ytitle = "d#sigma(WZ#rightarrow3l#nu)/d" + variable + " (pb)";
    }
  else
    {
      xtitle = "p_{T}^{" + variable + "} (GeV)";
      ytitle = "d#sigma(WZ#rightarrow3l#nu)/dp_{T}^{" + variable + "} (pb/GeV)";
    }

  TFile* file = new TFile("rootfiles/all_unfolding_" + var + ".root", "read");

  TH1D* xsValue          = (TH1D*)file->Get("hComb_diff");
  TH1D* xsValue_Madgraph = (TH1D*)file->Get("hGenXs" + var + "_1");
  TH1D* xsValue_MCnlo;

  if (drawMcfm) xsValue_MCnlo = (TH1D*)file->Get("mcfm_tot");


  // Data cosmetics
  //----------------------------------------------------------------------------
  xsValue->SetFillColor  (     kWhite);
  xsValue->SetFillStyle  (       1001);
  xsValue->SetLineColor  (     kBlack);
  xsValue->SetLineWidth  (          1);
  xsValue->SetMarkerColor(     kBlack);
  xsValue->SetMarkerSize (     _msize);
  xsValue->SetMarkerStyle(kFullCircle);


  // Madgraph cosmetics
  //----------------------------------------------------------------------------
  xsValue_Madgraph->SetFillColor  (    kOrange);
  xsValue_Madgraph->SetFillStyle  (       1001);
  xsValue_Madgraph->SetLineColor  (  kOrange+7);
  xsValue_Madgraph->SetLineWidth  (          1);
  xsValue_Madgraph->SetMarkerColor(  kOrange+7);
  xsValue_Madgraph->SetMarkerSize (     _msize);
  xsValue_Madgraph->SetMarkerStyle(kFullSquare);


  // MCNLO cosmetics
  //----------------------------------------------------------------------------
  if (drawMcfm)
    {
      xsValue_MCnlo->SetFillColor  (   kAzure-9);
      xsValue_MCnlo->SetFillStyle  (       1001);
      xsValue_MCnlo->SetLineColor  (     kAzure);
      xsValue_MCnlo->SetLineWidth  (          1);
      xsValue_MCnlo->SetMarkerColor(     kAzure);
      xsValue_MCnlo->SetMarkerSize (     _msize);
      xsValue_MCnlo->SetMarkerStyle(kOpenCircle);
    }


  // Set the canvas and pads
  //----------------------------------------------------------------------------
  TCanvas* canvas = new TCanvas(var, var);

  TPad* pad1;
  TPad* pad2;
  TPad* pad3;

  if (drawMcfm)
    {
      pad1 = new TPad("pad1" + var, "pad1" + var, 0, 0.49, 1, 1.000);
      pad2 = new TPad("pad2" + var, "pad2" + var, 0, 0.33, 1, 0.492);
      pad3 = new TPad("pad3" + var, "pad3" + var, 0, 0.00, 1, 0.332);
  
      pad1->SetTopMargin(0.09);
      pad2->SetTopMargin(0);
      pad3->SetTopMargin(0);

//.........这里部分代码省略.........

//.........这里部分代码省略.........
	TH1D * vhCorrPtRef[2][5], *vhCorrPt[2][5];
	Int_t colors[10] = {kBlack,kRed,kYellow+2,kGreen+2,kBlue};
   Int_t styles[2] = {kFullCircle,kOpenCircle};

   int icent=0;
 //  int icentRef=numCentBinRef-1;

	for (Int_t lv=0; lv<2; ++lv) {
		for (Int_t c=0; c<numCentBin; ++c) {
			vhCorrPt[lv][c] = (TH1D*) trkCorr->InspectCorr(lv,c,c,jetBegBin,jetEndBin,2,7-etaPM-1,7+etaPM);
			handsomeTH1(vhCorrPt[lv][c],kBlue,1,kOpenCircle);
         vhCorrPt[lv][icent]->SetAxisRange(xmin,xmax,"X");
         vhCorrPt[lv][icent]->SetAxisRange(0,1,"Y");
		}
		for (Int_t c=0; c<numCentBinRef; ++c) {
			vhCorrPtRef[lv][c] = (TH1D*) trkCorrRef->InspectCorr(lv,c,c,jetBegBin,jetEndBin,2,7-etaPM-1,7+etaPM);
			handsomeTH1(vhCorrPtRef[lv][c],colors[c]);
         vhCorrPtRef[lv][icent]->SetAxisRange(xmin,xmax,"X");
         vhCorrPtRef[lv][icent]->SetAxisRange(0,1,"Y");
		}
	}
   
	TCanvas * cEff = new TCanvas("cEff","cEff",500,500);
   cEff->SetLogx();
   vhCorrPt[0][icent]->SetTitle(";Track p_{T} (GeV/c);A #times #epsilon");
   vhCorrPt[0][icent]->SetTitleOffset(1.2);
   vhCorrPt[0][icent]->SetTitleSize(0.055);
	vhCorrPt[0][icent]->Draw("E");
	vhCorrPt[1][icent]->Draw("sameE");
	vhCorrPtRef[0][0]->Draw("sameE");
	vhCorrPtRef[1][0]->Draw("sameE");
                for (Int_t c=0; c<numCentBinRef; ++c) {
	vhCorrPtRef[0][c]->Draw("sameE");
	vhCorrPtRef[1][c]->Draw("sameE");
}
   TLegend *leg0 = new TLegend(0.16,0.84,0.46,0.92);
   leg0->SetFillStyle(0);
   leg0->SetBorderSize(0);
   leg0->SetTextSize(0.04);
   leg0->AddEntry(vhCorrPt[0][0],"PYTHIA+HYDJET","");
   if (etaPM==5)   leg0->AddEntry(vhCorrPt[0][0],Form("Track |#eta|<2.4"),"");
   if (etaPM==2) leg0->AddEntry(vhCorrPt[0][0],Form("Track |#eta|<1.2"),"");
	leg0->Draw();
   TLine * l = new TLine(xmin,1,xmax,1);
   l->SetLineStyle(2);
   l->Draw();
	
   TLegend *leg = new TLegend(0.34,0.32,0.60,0.48);
   leg->SetFillStyle(0);
   leg->SetBorderSize(0);
   leg->SetTextSize(0.035);
   leg->AddEntry(vhCorrPt[0][icent],title,"p");
    if (numCentBinRef==2) {
        leg->AddEntry(vhCorrPtRef[0][0],reftitle+ "0-30%","p");
        leg->AddEntry(vhCorrPtRef[0][1],reftitle+ "30-100%","p");
    } else {
        leg->AddEntry(vhCorrPtRef[0][0],reftitle+ "0-10%","p");
        leg->AddEntry(vhCorrPtRef[0][1],reftitle+ "10-30%","p");
        leg->AddEntry(vhCorrPtRef[0][2],reftitle+ "30-50%","p");
        leg->AddEntry(vhCorrPtRef[0][3],reftitle+ "50-100%","p");
    }
   leg->Draw();
   
	drawText("CMS Simulation",0.64,0.89);
	drawText("Fake Rate",0.69,0.26);
   
   cEff->Print(outdir+"/"+tag+".gif");
   cEff->Print(outdir+"/"+tag+".pdf");

	TCanvas * cJet = new TCanvas("cJet","cJet",500,500);
	cJet->SetLogy();

	trkCorr->vhPtHat[1][icent]->Draw("E");
        
       for (Int_t c=0; c<numCentBinRef; ++c) {
	trkCorrRef->vhPtHat[1][c]->SetMarkerColor(kRed+c);
	trkCorrRef->vhPtHat[1][c]->SetMarkerStyle(kOpenCircle+c);
	trkCorrRef->vhPtHat[1][c]->Draw("same E");
        }

	TCanvas * cCent = new TCanvas("cCent","cCent",500,500);
   TH1D * hCent = (TH1D*)trkCorr->sample_[0]->Get("hCent");
   TH1D * hCentRef = (TH1D*)trkCorrRef->sample_[0]->Get("hCent");
	hCentRef->SetMarkerStyle(kOpenCircle);
	hCentRef->Scale(1./hCentRef->GetEntries());
	hCent->Scale(1./hCent->GetEntries());
	hCent->Draw("p");
	hCentRef->Draw("same p");

    if(SaveFile){
        TFile * outf = new TFile(Form("%s/%sTrkEffFake.root", outdir.Data(),tag.Data()), "CREATE");
        for (Int_t lv=0; lv<2; ++lv) {
        vhCorrPt[lv][icent]->Write();
            for (Int_t c=numCentBinRef-1; c>=0; --c) {
                vhCorrPtRef[lv][c]->Write();
            }
        }
        outf->Close();
    }
}

void PtSpectraCombination_StagedNormalization_EtaCM_M08_M13_NoOLDAlignmentRuns()
{
   gROOT->Reset();
   gROOT->ProcessLine(".x rootlogonChristof.C");
   gROOT->ForceStyle();
   gStyle->SetPalette(1);
   gStyle->SetTitleYOffset(1.27);
   gStyle->SetOptFit(0000);
//   gStyle->SetPadRightMargin(0.16);

   bool doSave = true;

   TFile *f_minbias = new TFile("CombineSpectra_minbias_EtaCM_M08_M13_NoOLDAlignmentRuns.root");
   TFile *f_FullTrack_12 = new TFile("CombineSpectra_FullTrackTrigger_Track12_EtaCM_M08_M13_NoOLDAlignmentRuns.root");
   TFile *f_FullTrack_20 = new TFile("CombineSpectra_FullTrackTrigger_Track20_EtaCM_M08_M13_NoOLDAlignmentRuns.root");
   TFile *f_FullTrack_30 = new TFile("CombineSpectra_FullTrackTrigger_Track30_EtaCM_M08_M13_NoOLDAlignmentRuns.root");
   TFile *f_out = new TFile("PtSpectraCombination_StagedNormalization_EtaCM_M08_M13_NoOLDAlignmentRuns_TrackTrigger.root","recreate");

   //Get spectra
   TH1D* hPartPt_minbias_corrected = (TH1D*)f_minbias->Get("hPartPt_minbias_trkCorr_trigCorr");
   TH1D* hPartPt_0_14_minbias_corrected = (TH1D*)f_minbias->Get("hPartPt_0_14_minbias_trkCorr_trigCorr");
   TH1D* hPartPt_14_22_FullTrack12_corrected = (TH1D*)f_FullTrack_12->Get("hPartPt_FullTrack12_14_22_trkCorr_trigCorr");
   TH1D* hPartPt_22_32_FullTrack20_corrected = (TH1D*)f_FullTrack_20->Get("hPartPt_FullTrack20_22_32_trkCorr_trigCorr");
   TH1D* hPartPt_32_X_FullTrack30_corrected = (TH1D*)f_FullTrack_30->Get("hPartPt_FullTrack30_SpectComb_trkCorr_trigCorr");

   hPartPt_0_14_minbias_corrected->SetBinContent(hPartPt_0_14_minbias_corrected->FindBin(103.24),0);//Remove muon event
   hPartPt_0_14_minbias_corrected->SetBinError(hPartPt_0_14_minbias_corrected->FindBin(103.24),0);//Remove muon event
   hPartPt_minbias_corrected->SetBinContent(hPartPt_minbias_corrected->FindBin(103.24),0);//Remove muon event
   hPartPt_minbias_corrected->SetBinError(hPartPt_minbias_corrected->FindBin(103.24),0);//Remove muon event

   //Get number of events
   TH1D* hNumEv_minbias = (TH1D*)f_minbias->Get("hNumEv_minbias_trigCorr");
   TH1D* hNumEv_0_14_minbias = (TH1D*)f_minbias->Get("hNumEv_0_14_minbias_trigCorr");
   TH1D* hNumEv_14_22_minbias = (TH1D*)f_minbias->Get("hNumEv_14_22_minbias_trigCorr");
   TH1D* hNumEv_14_22_FullTrack12 = (TH1D*)f_FullTrack_12->Get("hNumEv_FullTrack12_14_22_trigCorr");
   TH1D* hNumEv_22_32_FullTrack12 = (TH1D*)f_FullTrack_12->Get("hNumEv_FullTrack12_22_32_trigCorr");
   TH1D* hNumEv_22_32_FullTrack20 = (TH1D*)f_FullTrack_20->Get("hNumEv_FullTrack20_22_32_trigCorr");
   TH1D* hNumEv_32_X_FullTrack20 = (TH1D*)f_FullTrack_20->Get("hNumEv_FullTrack20_32_X_trigCorr");
   TH1D* hNumEv_32_X_FullTrack30 = (TH1D*)f_FullTrack_30->Get("hNumEv_FullTrack30_32_X_trigCorr");

   float numev_minbias = hNumEv_minbias->GetBinContent(1);
   float numev_0_14_minbias = hNumEv_0_14_minbias->GetBinContent(1);
   float numev_14_22_minbias = hNumEv_14_22_minbias->GetBinContent(1);
   float numev_14_22_FullTrack12 = hNumEv_14_22_FullTrack12->GetBinContent(1);
   float numev_22_32_FullTrack12 = hNumEv_22_32_FullTrack12->GetBinContent(1);
   float numev_22_32_FullTrack20 = hNumEv_22_32_FullTrack20->GetBinContent(1);
   float numev_32_X_FullTrack20 = hNumEv_32_X_FullTrack20->GetBinContent(1);
   float numev_32_X_FullTrack30 = hNumEv_32_X_FullTrack30->GetBinContent(1);

   std::cerr<<" Number of minbias events in the normalization classes:" << std::endl;
   std::cerr<<"   All      : " << numev_minbias << std::endl;
   std::cerr<<"    0-14 MB : " << numev_0_14_minbias << std::endl;
   std::cerr<<"   14-22 MB : " << numev_14_22_minbias << std::endl;
   std::cerr<<"   14-22 T12: " << numev_14_22_FullTrack12 << std::endl;
   std::cerr<<"   22-32 T12: " << numev_22_32_FullTrack12 << std::endl;
   std::cerr<<"   22-32 T20: " << numev_22_32_FullTrack20 << std::endl;
   std::cerr<<"   32-X  T20: " << numev_32_X_FullTrack20 << std::endl;
   std::cerr<<"   32-X  T30: " << numev_32_X_FullTrack30 << std::endl;

   //Prepare the individual normalizations
   float norm_0_14_minbias = numev_minbias;
   float norm_14_22_FullTrack12 = numev_14_22_FullTrack12*numev_minbias/numev_14_22_minbias;
//   float norm_14_22_FullTrack12 = numev_14_22_FullTrack12*11165.040;
   float norm_22_32_FullTrack20 = norm_14_22_FullTrack12*numev_22_32_FullTrack20/numev_22_32_FullTrack12;
   float norm_32_X_FullTrack30 = norm_22_32_FullTrack20*numev_32_X_FullTrack30/numev_32_X_FullTrack20;

   std::cerr<<"Normalization factors:" << std::endl;
   std::cerr<<"      norm_0_14_minbias     : " << norm_0_14_minbias << std::endl;
   std::cerr<<"      norm_14_22_FullTrack12: " << norm_14_22_FullTrack12 << std::endl;
   std::cerr<<"      norm_22_32_FullTrack20: " << norm_22_32_FullTrack20 << std::endl;
   std::cerr<<"      norm_32_X_FullTrack30 : " << norm_32_X_FullTrack30 << std::endl;

   hPartPt_minbias_corrected->Scale(1./numev_minbias);
   hPartPt_0_14_minbias_corrected->Scale(1./norm_0_14_minbias);
   hPartPt_14_22_FullTrack12_corrected->Scale(1./norm_14_22_FullTrack12);
   hPartPt_22_32_FullTrack20_corrected->Scale(1./norm_22_32_FullTrack20);
   hPartPt_32_X_FullTrack30_corrected->Scale(1./norm_32_X_FullTrack30);

   TCanvas *c1 = new TCanvas("c1","c1");
   c1->cd();
   c1->SetLogy();
   hPartPt_0_14_minbias_corrected->GetXaxis()->CenterTitle();
   hPartPt_0_14_minbias_corrected->GetYaxis()->CenterTitle();
   hPartPt_0_14_minbias_corrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
   hPartPt_0_14_minbias_corrected->GetYaxis()->SetTitle("1/N_{ev} dN/dp_{T} (|#eta_{CM}|<1)");
//   hPartPt_0_14_minbias_corrected->SetNdivisions(505);
   hPartPt_0_14_minbias_corrected->SetMarkerColor(42);
   hPartPt_0_14_minbias_corrected->SetLineColor(42);
   hPartPt_0_14_minbias_corrected->SetMarkerStyle(34);
   hPartPt_0_14_minbias_corrected->SetMinimum(1e-11);
   hPartPt_0_14_minbias_corrected->SetFillColor(42);
//   hPartPt_0_14_minbias_corrected->SetFillStyle(3305);
   hPartPt_0_14_minbias_corrected->Draw("hist");

   hPartPt_minbias_corrected->SetMarkerStyle(20);
   hPartPt_minbias_corrected->Draw("same");

   hPartPt_14_22_FullTrack12_corrected->SetMarkerColor(2);
   hPartPt_14_22_FullTrack12_corrected->SetLineColor(2);
   hPartPt_14_22_FullTrack12_corrected->SetMarkerStyle(20);
//.........这里部分代码省略.........

TH1D* getYield(TTree* nt, TTree* ntMC, TString triggerpass, TString triggername, TString prescale, TString variable, TString varname, TString varlatex, Int_t BIN_NUM, Double_t BIN_MIN, Double_t BIN_MAX, TString addcut="")
{
  TH1D* hDistrib = new TH1D(Form("h%s_Distrib_%s",triggername.Data(),varname.Data()),"",BIN_NUM,BIN_MIN,BIN_MAX);
  for(float ivar=0;ivar<BIN_NUM;ivar++)
  {
    TCanvas* c = new TCanvas(Form("c%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"",500,500);
    TH1D* h = new TH1D(Form("h%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),";D^{0} mass (GeV/c^{2});Candidates",60,1.7,2.0);
    TH1D* hMC = new TH1D(Form("hMC%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"",60,1.75,1.95);
    TH1D* hSW = new TH1D(Form("hSW%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"",60,1.75,1.95);
    Float_t varmin = BIN_MIN+ivar*((BIN_MAX-BIN_MIN)/BIN_NUM);
    Float_t varmax = BIN_MIN+(ivar+1)*((BIN_MAX-BIN_MIN)/BIN_NUM);
    nt->Project(Form("h%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),Form("Dmass%s",prescale.Data()),Form("%s%s&&(%s>%f&&%s<%f)%s",prefilter.Data(),addcut.Data(),variable.Data(),varmin,variable.Data(),varmax,triggerpass.Data()));
    ntMC->Project(Form("hMC%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"Dmass",Form("%s%s&&(%s>%f&&%s<%f)%s",prefilterMC.Data(),addcut.Data(),variable.Data(),varmin,variable.Data(),varmax,"&&1"));
    ntMC->Project(Form("hSW%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"Dmass",Form("%s%s&&(%s>%f&&%s<%f)%s",prefilterSW.Data(),addcut.Data(),variable.Data(),varmin,variable.Data(),varmax,"&&1"));
    h->SetMaximum(h->GetMaximum()*1.20);
    h->Draw();

    TF1* f = new TF1(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"[0]*([7]*([9]*Gaus(x,[1],[2]*(1+[11]))/(sqrt(2*3.14159)*[2]*(1+[11]))+(1-[9])*Gaus(x,[1],[10]*(1+[11]))/(sqrt(2*3.14159)*[10]*(1+[11])))+(1-[7])*Gaus(x,[1],[8]*(1+[11]))/(sqrt(2*3.14159)*[8]*(1+[11])))+[3]+[4]*x+[5]*x*x+[6]*x*x*x", 1.7, 2.0);
    f->SetParLimits(4,-1000,1000);
    f->SetParLimits(10,0.001,0.05);
    f->SetParLimits(2,0.01,0.1);
    f->SetParLimits(8,0.02,0.2);
    f->SetParLimits(7,0,1);
    f->SetParLimits(9,0,1);

    f->SetParameter(0,setparam0);
    f->SetParameter(1,setparam1);
    f->SetParameter(2,setparam2);
    f->SetParameter(10,setparam10);
    f->SetParameter(9,setparam9);

    f->FixParameter(8,setparam8);
    f->FixParameter(7,1);
    f->FixParameter(1,fixparam1);
    f->FixParameter(3,0);
    f->FixParameter(4,0);
    f->FixParameter(5,0);
    f->FixParameter(6,0);
    f->FixParameter(11,0);
    h->GetEntries();

    hMC->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"q","",1.7,2.0);
    hMC->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"q","",1.7,2.0);
    f->ReleaseParameter(1);
    hMC->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
    hMC->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
    hMC->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L m","",1.7,2.0);

    f->FixParameter(1,f->GetParameter(1));
    f->FixParameter(2,f->GetParameter(2));
    f->FixParameter(10,f->GetParameter(10));
    f->FixParameter(9,f->GetParameter(9));
    f->FixParameter(7,0);
    f->ReleaseParameter(8);
    f->SetParameter(8,setparam8);

    hSW->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
    hSW->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
    hSW->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
    hSW->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L m","",1.7,2.0);

    f->FixParameter(7,hMC->Integral(0,1000)/(hSW->Integral(0,1000)+hMC->Integral(0,1000)));
    f->FixParameter(8,f->GetParameter(8));
    f->ReleaseParameter(3);
    f->ReleaseParameter(4);
    f->ReleaseParameter(5);
    f->ReleaseParameter(6);

    f->SetLineColor(kRed);

    h->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"q","",1.7,2.0);
    h->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"q","",1.7,2.0);
    f->ReleaseParameter(1);
    f->SetParLimits(1,1.86,1.87);
    f->ReleaseParameter(11);
    f->SetParLimits(11,-1.,1.);
    h->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
    h->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
    h->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
    h->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L m","",1.7,2.0);
    h->SetMarkerSize(0.8);
    h->SetMarkerStyle(20);

    TF1* background = new TF1(Form("background%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"[0]+[1]*x+[2]*x*x+[3]*x*x*x");
    background->SetParameter(0,f->GetParameter(3));
    background->SetParameter(1,f->GetParameter(4));
    background->SetParameter(2,f->GetParameter(5));
    background->SetParameter(3,f->GetParameter(6));
    background->SetLineColor(4);
    background->SetRange(1.7,2.0);
    background->SetLineStyle(2);

    TF1* mass = new TF1(Form("fmass%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"[0]*([3]*([4]*Gaus(x,[1],[2]*(1+[6]))/(sqrt(2*3.14159)*[2]*(1+[6]))+(1-[4])*Gaus(x,[1],[5]*(1+[6]))/(sqrt(2*3.14159)*[5]*(1+[6]))))");
    mass->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2),f->GetParameter(7),f->GetParameter(9),f->GetParameter(10),f->GetParameter(11));
    mass->SetParError(0,f->GetParError(0));
    mass->SetParError(1,f->GetParError(1));
    mass->SetParError(2,f->GetParError(2));
    mass->SetParError(3,f->GetParError(7));
    mass->SetParError(4,f->GetParError(9));
    mass->SetParError(5,f->GetParError(10));
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  RooDataHist ewkMet ("ewkMET", "ewkMET", RooArgSet(pfmet),hEWKMet);  RooHistPdf pdfEWK ("ewk", "ewk", pfmet,ewkMet, 1);
  RooDataHist ewkMetp("ewkMETp","ewkMETp",RooArgSet(pfmet),hEWKMetp); RooHistPdf pdfEWKp("ewkp","ewkp",pfmet,ewkMetp,1); 
  RooDataHist ewkMetm("ewkMETm","ewkMETm",RooArgSet(pfmet),hEWKMetm); RooHistPdf pdfEWKm("ewkm","ewkm",pfmet,ewkMetm,1); 
  
  // QCD Pdfs
  CPepeModel1 qcd("qcd",pfmet);
  CPepeModel1 qcdp("qcdp",pfmet);
  CPepeModel1 qcdm("qcdm",pfmet);
  
  // Signal + Background PDFs
  RooAddPdf pdfMet ("pdfMet", "pdfMet", RooArgList(pdfWe,pdfEWK,*(qcd.model)),   RooArgList(nSig,nEWK,nQCD));  
  RooAddPdf pdfMetp("pdfMetp","pdfMetp",RooArgList(pdfWep,pdfEWKp,*(qcdp.model)),RooArgList(nSigp,nEWKp,nQCDp));
  RooAddPdf pdfMetm("pdfMetm","pdfMetm",RooArgList(pdfWem,pdfEWKm,*(qcdm.model)),RooArgList(nSigm,nEWKm,nQCDm));
    
  //
  // Perform fits
  //
  RooDataHist dataMet("dataMet", "dataMet", RooArgSet(pfmet),hDataMet);
  RooFitResult *fitRes = pdfMet.fitTo(dataMet,Extended(),Minos(kTRUE),Save(kTRUE));
  
  RooDataHist dataMetp("dataMetp","dataMetp",RooArgSet(pfmet),hDataMetp);
  RooFitResult *fitResp = pdfMetp.fitTo(dataMetp,Extended(),Minos(kTRUE),Save(kTRUE));
  
  RooDataHist dataMetm("dataMetm","dataMetm",RooArgSet(pfmet),hDataMetm);
  RooFitResult *fitResm = pdfMetm.fitTo(dataMetm,Extended(),Minos(kTRUE),Save(kTRUE));
    
  //
  // Use histogram version of fitted PDFs to make ratio plots
  // (Will also use PDF histograms later for Chi^2 and KS tests)
  //
  TH1D *hPdfMet = (TH1D*)(pdfMet.createHistogram("hPdfMet", pfmet));
  hPdfMet->Scale((nSig.getVal()+nEWK.getVal()+nQCD.getVal())/hPdfMet->Integral());
  TH1D *hMetDiff = makeDiffHist(hDataMet,hPdfMet,"hMetDiff");
  hMetDiff->SetMarkerStyle(kFullCircle);
  hMetDiff->SetMarkerSize(0.9);
   
  TH1D *hPdfMetp = (TH1D*)(pdfMetp.createHistogram("hPdfMetp", pfmet));
  hPdfMetp->Scale((nSigp.getVal()+nEWKp.getVal()+nQCDp.getVal())/hPdfMetp->Integral());
  TH1D *hMetpDiff = makeDiffHist(hDataMetp,hPdfMetp,"hMetpDiff");
  hMetpDiff->SetMarkerStyle(kFullCircle);
  hMetpDiff->SetMarkerSize(0.9);
    
  TH1D *hPdfMetm = (TH1D*)(pdfMetm.createHistogram("hPdfMetm", pfmet));
  hPdfMetm->Scale((nSigm.getVal()+nEWKm.getVal()+nQCDm.getVal())/hPdfMetm->Integral());
  TH1D *hMetmDiff = makeDiffHist(hDataMetm,hPdfMetm,"hMetmDiff");
  hMetmDiff->SetMarkerStyle(kFullCircle); 
  hMetmDiff->SetMarkerSize(0.9);
   
  
  //--------------------------------------------------------------------------------------------------------------
  // Make plots 
  //==============================================================================================================  
  
  TCanvas *c = MakeCanvas("c","c",800,800);
  c->Divide(1,2,0,0);
  c->cd(1)->SetPad(0,0.3,1.0,1.0);
  c->cd(1)->SetTopMargin(0.1);
  c->cd(1)->SetBottomMargin(0.01);
  c->cd(1)->SetLeftMargin(0.18);  
  c->cd(1)->SetRightMargin(0.07);  
  c->cd(1)->SetTickx(1);
  c->cd(1)->SetTicky(1);  
  c->cd(2)->SetPad(0,0,1.0,0.3);
  c->cd(2)->SetTopMargin(0.05);
  c->cd(2)->SetBottomMargin(0.45);
  c->cd(2)->SetLeftMargin(0.18);

//.........这里部分代码省略.........
  hs->Add(qcd1);
  hs->Add(qcd2);
  hs->Add(qcd3);
  hs->Add(qcd4);
  hs->Add(qcd5);
  hs->Add(qcd6);
  hs->Add(qcd7);
  hs->Add(qcd8);
  hs->Add(qcd9);
  hs->Add(qcd10);
  hs->Add(qcd11);
  }
  
  
  if(inclZ == true){
  hs->Add(zjets);
  }else{
  hs->Add(z1jets);
  hs->Add(z2jets);
  hs->Add(z3jets);
  hs->Add(z4jets);  
  }
  
  if(inclW == true){
  hs->Add(wjets);
  }else{
  hs->Add(w1jets);
  hs->Add(w2jets);
  hs->Add(w3jets);
  hs->Add(w4jets);  
  }
      
  hs->Add(top_t);
  hs->Add(top_tw);
  hs->Add(top_s);
  hs->Add(tbar_t);
  hs->Add(tbar_tw);
  hs->Add(tbar_s);
  
  hs->Add(tt);

  //draw histos to files
  TCanvas *c1 = new TCanvas("Plot","Plot",900, 600);
		
  hs->SetMaximum(data->GetBinContent(data->GetMaximumBin())*1.2);

  hs->Draw();
  data->Draw("E same");
  data->SetMarkerStyle(20);

//events:
cout << "ttbar: " << tt->Integral() << endl;
cout << "data: " << data->Integral() << endl;
  
  hs->GetXaxis()->SetLimits(MinX, MaxX);
  hs->GetXaxis()->SetTitle(Xtitle); hs->GetXaxis()->SetTitleSize(0.05);
  hs->GetYaxis()->SetTitle("Number of Events");hs->GetYaxis()->SetTitleSize(0.05);
  
  
  	TLegend *tleg2;
	tleg2 = new TLegend(0.7,0.7,0.8,0.9);
	tleg2->SetTextSize(0.04);
	tleg2->SetBorderSize(0);
	tleg2->SetFillColor(10);
	tleg2->AddEntry(data , "2012 data", "lpe");
	tleg2->AddEntry(tt , "t#bar{t}", "lf");
	tleg2->AddEntry(top_t, "single top", "lf");
	tleg2->AddEntry(wjets , "w+jets", "lf");
	tleg2->AddEntry(zjets , "z+jets", "lf");
	tleg2->AddEntry(qcd , "QCD", "lf");
	
	//tleg2->AddEntry(singtEff, "single-t"      , "l");
	//tleg2->AddEntry(singtwEff, "single-tW"      , "l");
 	tleg2->Draw("same");	
	
	TText* textPrelim = doPrelim(0.17,0.96);
	textPrelim->Draw();
	
  if(logPlot ==true){
  c1->SetLogy();
  }	
  
  TString plotName("plots/Control/Btags/");
  
  if(logPlot ==true){
    plotName += Variable+"Test_Log";
    plotName += ".pdf";
    
  }else{
    plotName += Variable;  
    plotName += ".pdf";
  }
 
 
  c1->SaveAs(plotName);
  delete c1;
  
  }
  	
}

//.........这里部分代码省略.........
   
  ntMix->Draw("var>>hDataMixFull",weight*(dijet&&jetID&&centHF&&vtx));
  ntMix->Draw("var>>hDataMixB",weight*(side&&jetID&&centHF&&vtx));
  ntMix->Draw("pt1>>hNormDataMix",weight*(lead&&jetID&&centHF&&vtx));
  
  ntReference->Draw("var>>hReferenceFull",dijet&&noise&&jetID&&centHF&&vtx);
  ntReference->Draw("var>>hReferenceB",side&&noise&&jetID&&centHF&&vtx);
  ntReference->Draw("pt1>>hNormReference",lead&&noise&&jetID&&centHF&&vtx);

  hDataMix->Add(hDataMixFull);
  h->Add(hFull);
  hReference->Add(hReferenceFull);

  if(subtract){
    hDataMix->Add(hDataMixB,-1);
    h->Add(hB,-1);
    hReference->Add(hReferenceB,-1);
  }

  hB->SetFillStyle(3005);
  hB->SetFillColor(15);

  // calculate the statistical error and normalize
  h->SetLineColor(dataColor);
  h->SetMarkerColor(dataColor);
  h->Sumw2();
  if(normLead){
    h->Scale(1./hNorm->Integral());
    hB->Scale(1./hNorm->Integral());
  }else{
    hB->Scale(1./h->Integral());
    h->Scale(1./h->Integral());
  }
  h->SetMarkerStyle(20);

  if(hReference->Integral() > 0){
    hReference->Scale(1./hNormReference->Integral());
  }

  hReference->SetLineColor(kBlue);
  hReference->SetFillColor(kAzure-8);
  hReference->SetFillStyle(3005);

  if(normLead){
    hDataMixB->Scale(1./hNormDataMix->Integral());
    hDataMix->Scale(1./hNormDataMix->Integral());
  }else{
    hDataMixB->Scale(1./hDataMix->Integral());
    hDataMix->Scale(1./hDataMix->Integral());
  }
  hDataMix->SetLineColor(mixColor);
  hDataMix->SetFillColor(mixColor);
  hDataMix->SetFillStyle(3004);
  
  hDataMix->SetMarkerSize(0);
  hDataMix->SetStats(0);

  hDataMix->GetXaxis()->SetLabelSize(22);
  hDataMix->GetXaxis()->SetLabelFont(43);
  hDataMix->GetXaxis()->SetTitleSize(28);
  hDataMix->GetXaxis()->SetTitleFont(43);
  hDataMix->GetXaxis()->SetTitleOffset(2.2);
  hDataMix->GetXaxis()->CenterTitle();
  
  
  hDataMix->GetYaxis()->SetLabelSize(22);

void doControlPlotsMET(){
setTDRStyle();

//loop over variables
for(int i = 0; i<N; i++){
//double MinX = MinXs[i];
//double MaxX = MaxXs[i];
Variable = Variables[i];
TString Xtitle = XTitles[i];
int RebingFact = RebinFacts[i];

//Data
TH1D* data = getSample("DoubleMu", 1, Obj, Next, Variable, RebinFact, Systematic);

std::cout << data->Integral() << std::endl;

//MC
TH1D* tt = getSample("TTJet", 1, Obj, Next, Variable, RebinFact, Systematic);
TH1D* wjets = getSample("WJetsToLNu", 1, Obj, Next, Variable, RebinFact, Systematic);
TH1D* DY1 = getSample("DYJetsToLL_M-10To50", 1, Obj, Next, Variable, RebinFact, Systematic);
TH1D* DY2 = getSample("DYJetsToLL_M-50", 1, Obj, Next, Variable, RebinFact, Systematic);
TH1D* T_tW = getSample("T_tW-channel", 1, Obj, Next, Variable, RebinFact, Systematic);
TH1D* Tbar_tW = getSample("Tbar_tW-channel",1, Obj, Next, Variable, RebinFact, Systematic);

THStack *hs = new THStack("hs","test");

  hs->Add(wjets);
  hs->Add(DY1);
  hs->Add(DY2);
  hs->Add(T_tW); 
  hs->Add(Tbar_tW);
  hs->Add(tt); 
  
  //draw histos to files
  TCanvas *c1 = new TCanvas("Plot","Plot",900, 600);
		
  hs->SetMaximum(data->GetBinContent(data->GetMaximumBin())*1.3);

  hs->Draw();
  data->Draw("E same");
  data->SetMarkerStyle(20);
  
//  hs->GetXaxis()->SetLimits(MinX, MaxX);
  hs->GetXaxis()->SetTitle(Xtitle); hs->GetXaxis()->SetTitleSize(0.05);
  hs->GetYaxis()->SetTitle("Number of Events");hs->GetYaxis()->SetTitleSize(0.05);
   
   
  if(logPlot == true){
  c1->SetLogy();
  }	
  
  	TLegend *tleg2;
	tleg2 = new TLegend(0.6,0.7,0.8,0.9);
	tleg2->SetTextSize(0.04);
	tleg2->SetBorderSize(0);
	tleg2->SetFillColor(10);
	tleg2->AddEntry(data , "2012 data", "lpe");
	tleg2->AddEntry(tt , "t#bar{t}", "lf");
	tleg2->AddEntry(T_tW, "single-t", "lf");
	tleg2->AddEntry(Tbar_tW, "single-#bar{t}", "lf");	
	tleg2->AddEntry(DY1 , "DYJetsToLL-10To50", "lf");
	tleg2->AddEntry(DY2 , "DYJetsToLL-50", "lf");
	tleg2->AddEntry(wjets , "W+jets", "lf");
	
 	tleg2->Draw("same");	

  TString plotName("plots/Control/PassesCutsUpTo1Btag/MET/");
  
  if(logPlot == true){
    plotName += Variable+"Log.pdf";
//    plotName += Nbtags+".pdf";
    
  }else{
    plotName += Variable+".pdf";  
//    plotName += Nbtags+".pdf";
  }
 
  TText* textPrelim = doPrelim(0.16,0.96); 
  textPrelim->Draw();
  
  c1->SaveAs(plotName);
  delete c1;
  
  }
  	
}

void DoComparison(TString newFileName, TString oldFileName, vector<TString> newHistName, vector<TString> oldHistName)
{
  TFile newF(newFileName,"read");
  TDirectory *newDir = newF.GetDirectory("Merged");

  TFile old(oldFileName);

  TFile output("Comparison.root","update");

  assert(newHistName.size() == oldHistName.size());


  vector<TH1D*> newHists;
  vector<TH1D*> oldHists;
  
  for(UInt_t i = 0; i < newHistName.size(); i++) {
    TH1D *newH = (TH1D*)newDir->Get(newHistName[i]);
    if(!newH) {
      cout<<"Could not find "<<newHistName[i]<< " in "<<newDir->GetName()<<endl;
      return;
    }
    TH1D *oldH = (TH1D*)old.Get(oldHistName[i]);
    if(!oldH) {
      cout<<"Could not find "<<oldHistName[i]<< " in "<<old.GetName()<<endl;
      return;
    }
    
    newHists.push_back(newH);
    oldHists.push_back(oldH);
  }

  for(UInt_t i = 0; i < newHists.size(); i++) {

    // Make ratio
    TH1D *ratio = (TH1D*)newHists[i]->Clone();
    ratio->Divide(oldHists[i]);

    ratio->SetDirectory(0);
    TString newName = newHists[i]->GetName();
    newName += "Ratio";
    ratio->SetName(newName);
    ratio->SetTitle(newName);
    ratio->SetMarkerStyle(20);
    output.cd();
    ratio->Write(newName, TObject::kOverwrite);
    cout<<"Wrote "<<newName<<" to "<< output.GetName()<<endl;

    // Make difference
    TH1D *difference = (TH1D*)newHists[i]->Clone();
    difference->Add(oldHists[i],-1.);
    difference->SetDirectory(0);
    newName = newHists[i]->GetName();
    newName += "Difference";
    difference->SetName(newName);
    difference->SetTitle(newName);
    difference->SetMarkerStyle(20);
    output.cd();
    difference->Write(newName, TObject::kOverwrite);
    cout<<"Wrote "<<newName<<" to "<< output.GetName()<<endl;

    // Make side by side plot
    TString canName = "Canvas";
    canName += newHists[i]->GetName();
    canName += "Comparison";
    TCanvas *can = new TCanvas(canName, canName);
    newHists[i]->SetLineColor(kBlack);
    newHists[i]->SetMarkerStyle(20);
    newHists[i]->SetMarkerColor(kBlack);
    oldHists[i]->SetLineColor(kRed);
    oldHists[i]->SetMarkerStyle(20);
    oldHists[i]->SetMarkerColor(kRed);
    oldHists[i]->Draw();
    newHists[i]->Draw("same");
    can->Write(canName, TObject::kOverwrite);
    delete can; can = NULL;
  }

}

//.........这里部分代码省略.........
//  c->SetPad(0,0,1.0,1.0);
//  c->SetTopMargin(0.1);
//  c->SetBottomMargin(0.15);
//  c->SetLeftMargin(0.15);  
//  c->SetRightMargin(0.07);  
//  c->SetTickx(1);
//  c->SetTicky(1);  
  gStyle->SetTitleOffset(1.400,"Y");
  TGaxis::SetMaxDigits(3);
  char ylabel[100];  // string buffer for y-axis label
  // label for lumi
  char lumitext[100];
  if(lumi<0.1) sprintf(lumitext,"%.1f pb^{-1} (8 TeV)",lumi*1000.);
  else         sprintf(lumitext,"%.2f fb^{-1}  at  #sqrt{s} = 8 TeV",lumi);
  char CMStext[100];
  sprintf(CMStext,"#font[61]{CMS}");
  char Preliminarytext[100];
  sprintf(Preliminarytext,"#font[52]{Preliminary}");
  char binlabel[50];
  
  // plot colors
  Int_t linecolorW   = kOrange-3;
  Int_t fillcolorW   = kOrange-2;
  Int_t linecolorEWK = kOrange+10;
  Int_t fillcolorEWK = kOrange+7;
  Int_t linecolorQCD = kViolet+2;
  Int_t fillcolorQCD = kViolet-5;

  //
  // Dummy histograms for TLegend
  // (I can't figure out how to properly pass RooFit objects...)
  //
  TH1D *hDummyData = new TH1D("hDummyData","",0,0,10);
  hDummyData->SetMarkerStyle(kFullCircle);
  hDummyData->SetMarkerSize(0.9);
  
  TH1D *hDummyW = new TH1D("hDummyW","",0,0,10);
  hDummyW->SetLineColor(linecolorW);
  hDummyW->SetFillColor(fillcolorW);
  hDummyW->SetFillStyle(1001);
  
  TH1D *hDummyEWK = new TH1D("hDummyEWK","",0,0,10);
  hDummyEWK->SetLineColor(linecolorEWK);
  hDummyEWK->SetFillColor(fillcolorEWK);
  hDummyEWK->SetFillStyle(1001);
  
  TH1D *hDummyQCD = new TH1D("hDummyQCD","",0,0,10);
  hDummyQCD->SetLineColor(linecolorQCD);
  hDummyQCD->SetFillColor(fillcolorQCD);
  hDummyQCD->SetFillStyle(1001);

// Wpt distribution=========================
  CPlot* plotWPptLog;
  CPlot* plotWPptDiffLog;
  CPlot* plotWMptLog;
  CPlot* plotWMptDiffLog;

//W plus pt distribution
  TH1D* hWptMC_p = (TH1D*)hDYToTauTauP->Clone("hWptMC_p");
  hWptMC_p->Add(hTTJetsP);
  hWptMC_p->Add(hWToTauNuP);
  hWptMC_p->Add(hDYToMuMuP);
  hWptMC_p->Add(hQCDWPpt);
  hWptMC_p->Add(hSigWPpt);
 
  TH1D* hWPptDiffLog;

//_______________________________________________________________________
TH1D * GetPedestalHistoAndSystAndSubtractPed(Int_t binSystem, Int_t binAssoc,Int_t binMeson,TH1D *histo, TGraphAsymmErrors* gr,TGraphAsymmErrors *&grout, TString canvasname,TGraphAsymmErrors *&grbaseOut,TGraphAsymmErrors *&grv2Out){

  Double_t value = 0, pedestal=0;
  TGraphAsymmErrors *grBaseHelp,*grV2;
  grbaseOut=new TGraphAsymmErrors();
  grbaseOut->SetName(Form("grbaselineUncFull_%s_%s_%s",sets[binSystem].Data(),strmesonpt[binMeson].Data(),pthadron[binAssoc].Data()));

  Double_t xuncFull,errxuncFull;
  Double_t xuncv2,errxuncv2;
  Int_t bin,bingr;
  TString path = pedestalfilenames[binSystem][binAssoc];


  if(binSystem==0){//pp
    grV2=0x0;
    cout << "pp -->  Reading File from path: " << path << endl;
    
    TFile * file = TFile::Open(path.Data(),"READ");
    TCanvas* c=(TCanvas*)file->Get(canvasname.Data());
    TH1D* h = (TH1D*)c->GetListOfPrimitives()->FindObject("FinalTrendPedestal");
    grBaseHelp=(TGraphAsymmErrors*)c->GetListOfPrimitives()->FindObject("fBaselineVariationSystematicsPedestal");    

    if(isReflectedData){
      xuncFull=3.25;
      errxuncFull=0.075;
      xuncv2=-0.15; 
      errxuncv2=0.075;
    }
    else{
      xuncFull=4.83;
      errxuncFull=0.075;
      xuncv2=-1.7;
      errxuncv2=0.075;
    }

    bin=h->FindBin(mesonptcenter[binMeson]);
    bingr=GetBinGraph(mesonptcenter[binMeson],grBaseHelp);


    pedestal=h->GetBinContent(bin);
    Double_t x,y,erryl,erryh;
    Printf("histo: x=%f, graph: %f",h->GetBinCenter(bin),x);
    grBaseHelp->GetPoint(bingr,x,y);
    erryl=grBaseHelp->GetErrorYlow(bingr);
    erryh=grBaseHelp->GetErrorYhigh(bingr);
    
    grbaseOut->SetPoint(0,xuncFull,0);
    grbaseOut->SetPointError(0,errxuncFull,errxuncFull,erryl,erryh);
    if(grV2){
      grV2->GetPoint(bingr,x,y);
      erryl=grV2->GetErrorYlow(bingr);
      erryh=grV2->GetErrorYhigh(bingr);	
      grv2Out->SetPoint(0,xuncFull,0);
      grv2Out->SetPointError(0,errxuncFull,errxuncFull,erryl,erryh);
    }
  }

  cout<<"Baseline being subtracted for:"<<sets[binSystem].Data()<<strmesonpt[binMeson].Data()<<pthadron[binAssoc].Data()<<endl;
  
  grout=(TGraphAsymmErrors*)gr->Clone(Form("grSub_%s_%s_%s",sets[binSystem].Data(),strmesonpt[binMeson].Data(),pthadron[binAssoc].Data()));
  
  TString nameoutput = histo->GetName();
  nameoutput += "_subtr_";
  nameoutput += "pedestal";
  nameoutput += Form("_%s_%s_%s",sets[binSystem].Data(),strmesonpt[binMeson].Data(),pthadron[binAssoc].Data());
   
  TH1D * outputhisto = (TH1D*)histo->Clone(nameoutput.Data());
  outputhisto->Reset();
  outputhisto->SetStats(kFALSE);
    

  for(Int_t iBin = 1; iBin <= histo->GetNbinsX();iBin++){
        
        
    value = histo->GetBinContent(iBin);
    value -= pedestal;
      
    outputhisto->SetBinContent(iBin,value);
      
    outputhisto->SetBinError(iBin,histo->GetBinError(iBin));
    Double_t x,y,eyl,eyh;
    gr->GetPoint(iBin-1,x,y);
    eyl=gr->GetErrorYlow(iBin-1);
    eyh=gr->GetErrorYhigh(iBin-1);
    //cout<<x<<"  "<<y<<"  "<<eyl<<"  "<<eyh<<endl;
    grout->SetPoint(iBin-1,x,y-pedestal);

  }
  cout<<"sub -> "<<outputhisto->GetBinContent(5)<<endl;
  
  outputhisto->SetXTitle("#Delta#varphi (rad)");
  outputhisto->SetYTitle("#frac{1}{#it{N}_{D}} #frac{d#it{N}^{assoc}}{d#Delta#varphi} - baseline (rad^{-1})");
  outputhisto->GetYaxis()->CenterTitle();
  outputhisto->SetMarkerColor(kBlack);
  outputhisto->SetLineColor(kBlack);
  outputhisto->SetMarkerStyle(20);
  
  //  outputhisto->GetYaxis()->SetTitleOffset(1.5);
  //  outputhisto->GetYaxis()->SetTitleFont(42);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    TF1 *NBD_fun = new
    TF1("NBD_fun","[0]*TMath::Gamma(x+[1])/(TMath::Gamma(x+1)*TMath::Gamma([1]))*TMath::Power([2]/[1],x)/TMath::Power([2]/[1]+1,x+[1])",0,100);
    NBD_fun->SetParameter(0,1);	//[0]: Normalized constant
    NBD_fun->SetParameter(1,(*kbest)[0]);	//[1]: k value
    NBD_fun->SetParameter(2,(*mubest)[0]);	//[2]: mu value

    TTree *t = (TTree*) fGlauber->Get("nt_p_Pb");
    Long_t Nevent;

    Nevent = (Long_t) t->GetEntries();

    Long_t Ev;
    Int_t Bino;
    Double_t Para, Bi_Para, Mult;
    Float_t Ncoll;
    t->SetBranchAddress("Ncoll",&Ncoll);

    for(Ev=0; Ev<Nevent; Ev++) {
        if(Ev%100000==0)	 cout<<"Have run "<<Ev<<" events"<<endl;
        t->GetEntry(Ev);
        Para = 0; //make sure that Para doesn't accumulate through loops
        for(Bino=0; Bino<Ncoll; Bino++) {
            Bi_Para = NBD_fun->GetRandom();
            Para += Bi_Para;
        }
        histo_exp->Fill(Para);
    }
    Double_t SumEvent, scale;
    SumEvent = histo_exp->Integral(xbinmin,xbinmax);
    scale = 1/SumEvent;
    TH1D *histo_exp_norm = (TH1D*) histo_exp->Clone();
    histo_exp_norm->Scale(scale);

    TCanvas *c1 = new TCanvas();
    gStyle->SetOptStat(kFALSE);

    double hfbin[]= {0,1,2,3,4,6,8,10,13,16,20,25,30,40,55,70,90};
    int nhfbin = 16;
    rehisto_obs_norm = (TH1D*)histo_obs_norm->Rebin(nhfbin,"rehisto_obs_norm",hfbin);
    normalizeByBinWidth(rehisto_obs_norm);
    rehisto_exp_norm = (TH1D*)histo_exp_norm->Rebin(nhfbin,"rehisto_exp_norm",hfbin);
    normalizeByBinWidth(rehisto_exp_norm);
    TH1D* ratio = (TH1D*)rehisto_obs_norm->Clone("ratio");
    ratio->Divide(rehisto_exp_norm);
    ratio->SetMaximum(1.2);
    ratio->SetMinimum(0);
    ratio->GetXaxis()->SetTitle("HF #Sigma E_{T} |#eta|>4");
    ratio->GetYaxis()->SetTitle("ratio");

    TFile *fDSeff = TFile::Open("/afs/cern.ch/work/q/qixu/private/RpA/GlobalEvent/CentrDep/pPbHijing_EffCorr_forNBD.root");
    TFile *ftreff = TFile::Open("/afs/cern.ch/user/q/qixu/CMSSW_6_2_5/src/Centrality/Correction/pPbHist_Hijing_TrandEs.root");
    TH1D* hbef = (TH1D*)ftreff->Get("hHFEnergy4");
    TH1D* rehbef = (TH1D*)hbef->Rebin(nhfbin,"rehHFEnergy4",hfbin);
    TH1D* haft = (TH1D*)ftreff->Get("hHFEnergy4_tr");
    TH1D* rehaft = (TH1D*)haft->Rebin(nhfbin,"rehHFEnergy4_tr",hfbin);
    TGraphAsymmErrors *gtreff = new TGraphAsymmErrors();
    gtreff->BayesDivide(rehaft,rehbef);

    TGraphAsymmErrors *geff = (TGraphAsymmErrors*)fDSeff->Get("regEffHF4");
    for(int i=0; i<geff->GetN(); i++) {
        geff->SetPointEXlow(i,0);
        geff->SetPointEXhigh(i,0);
        gtreff->SetPointEXlow(i,0);
        gtreff->SetPointEXhigh(i,0);
    }

    ratio->SetTitle("");
    ratio->SetLineColor(1);
    ratio->SetMarkerStyle(24);
    ratio->SetMarkerColor(1);
    ratio->SetMarkerSize(1.5);
    ratio->Draw("P");

    geff->SetMarkerStyle(33);
    geff->SetMarkerColor(2);
    geff->SetMarkerSize(1.5);
    geff->Draw("Psame");

    gtreff->SetMarkerStyle(21);
    gtreff->SetMarkerColor(4);
    gtreff->SetMarkerSize(1.3);
    gtreff->Draw("Psame");

    TLegend *leg = new TLegend(0.60, 0.2, 0.78, 0.4);
    leg->SetFillColor(10);
    leg->SetFillStyle(0);
    leg->SetBorderSize(0.035);
    leg->SetTextFont(42);
    leg->SetTextSize(0.045);
    leg->AddEntry(ratio,"data/fit","p");
    leg->AddEntry(geff,"DS efficiency","p");
    leg->AddEntry(gtreff,"Event selection efficiency","p");
    leg->Draw("same");

    TLine *l = new TLine(0,1,90,1);
    l->SetLineStyle(2);
    l->Draw("same");
    c1->SaveAs(Form("%sratiovseff.png",dirname.Data()));

}